!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Originally created by James McDonagh at the University of Manchester 2015, in the Popelier group        !
! Components of the dihedral subroutine are acknowledged to Dr Tanja van Mourik University of St Andrews  !
! Licensed under MIT License                                        !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Originally created by James.
! Version 1.1
! CHANGE LOG
! Version 1  : Subroutines for calculating strutcural parameters
! Version 1.1: Modular format and incorporation in the Hermes program. Additional routines for jmol images
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

module module_subroutines

  use module_functions
  use module_constants

  implicit none
contains

  subroutine bondlength(atom_co_1,atom_co_2,bond_length)
    real, dimension(3), intent(in) :: atom_co_1, atom_co_2
    real, dimension(3) :: bond_vec
    real, intent(inout) :: bond_length

    bond_vec = bond_vector(atom_co_1,atom_co_2)
    bond_length = norm(bond_vec)

  end subroutine bondlength

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

  subroutine bondangle(atom_co_1,atom_co_2,atom_co_3,angle,angle_deg)

    real, dimension(3), intent(in) :: atom_co_1, atom_co_2, atom_co_3
    real, dimension(3) :: vec1, vec2
    real :: normvec1, normvec2, norms, dotvec1vec2, argument
    real,intent(inout) :: angle, angle_deg

    vec1 = bond_vector(atom_co_1, atom_co_2)
    vec2 = bond_vector(atom_co_3, atom_co_2)

    normvec1 = norm(vec1)
    normvec2 = norm(vec2)
    norms = normvec1 * normvec2

    dotvec1vec2 = vector_dot(vec1,vec2)
!    print*,'Dot product',dotvec1vec2               ! JMcD DIAGNOSTIC PRINT

    argument = dotvec1vec2 / norms
!   print*, argument                                ! JMcD DIAGNOSTIC PRINT

    angle = acos(argument)
    angle_deg = angle * radians_to_degrees

  end subroutine bondangle

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

  subroutine dihedral(atom_co_1,atom_co_2,atom_co_3,atom_co_4,dihedangle,dihedangle_deg)

    integer :: i
    real, dimension(3), intent(in) :: atom_co_1, atom_co_2, atom_co_3, atom_co_4
    real, dimension(3) :: b1, b2, b3, cross1, cross2, m
    real :: norm_b1, norm_b2, norm_b3, y, x
    real, intent(inout) :: dihedangle, dihedangle_deg
    Character (len=1) :: debug="F" ! Note debug is define on = T or off = F here

    If(debug=="T")then
       Print*, "In subroutine dihedral debug is turned on output will be verbose"
       Print*, "atom 1 coordinates are ", atom_co_1
       Print*, "atom 2 coordinates are ", atom_co_2
       Print*, "atom 3 coordinates are ", atom_co_3
       Print*, "atom 4 coordinates are ", atom_co_4
    end If

    b1=bond_vector(atom_co_1, atom_co_2)
    b2=bond_vector(atom_co_2, atom_co_3)
    b3=bond_vector(atom_co_3, atom_co_4)

    norm_b1=norm(b1)
    norm_b2=norm(b2)
    norm_b3=norm(b3)

    b1=normalise(b1,norm_b1)
    b2=normalise(b2,norm_b2)
    b3=normalise(b3,norm_b3)

    If(debug=="T")then
       Print*, "normalised bond vectors are b1 ", b1
       Print*, "normalised bond vectors are b2 ", b2
       Print*, "normalised bond vectors are b3 ", b3
    end If

    cross1=vector_cross(b1,b2)
    cross2=vector_cross(b2,b3)

    If(debug=="T")then
       Print*, "cross product vector 1 ", cross1
       Print*, "cross product vector 2 ", cross2
    end If

    m=vector_cross(cross1,b2)

    y=vector_dot(m,cross2)
    x=vector_dot(cross1,cross2)

    dihedangle=-atan2(y,x)
    dihedangle_deg=dihedangle*radians_to_degrees

  end subroutine dihedral

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

  subroutine StripSpaces(string)                      ! JMcD Acknowldge Jauch 27 November 2014 http ://stackoverflow.com/questions/27179549/removing-whitespace-in-string
    character(len=*) :: string
    integer :: stringLen
    integer :: last, actual

    stringLen = len (string)
    last = 1
    actual = 1

    do while (actual.lt.stringLen)
        if (string(last:last) == ' ') then
            actual = actual + 1
            string(last:last) = string(actual:actual)
            string(actual:actual) = ' '
        else
            last = last + 1
            if (actual.lt.last) &
                actual = last
        endif
    end do

  end subroutine StripSpaces

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

  subroutine sort_asend(array,toend,sorted)              ! JMcD Acknowledge Rossetta code.org Sorting algorithms/Bubble sort  http://rosettacode.org/wiki/Sorting_algorithms/Bubble_sort#Fortran
    ! Dummy arguments
    integer,intent(in) :: toend
    real, intent(in) :: array(1:toend)
    real, intent(inout) :: sorted(1:toend)
    ! Variables
    real :: temp
    logical :: swap
    integer :: i, j

    temp = 0
    sorted = 0
    i = 0
    j = 0
    sorted = array

    do j = toend, 1, -1
       swap = .FALSE.
       do i = 1, j
          if (sorted(i).gt.sorted(i + 1)) then
             temp = sorted(i)
             sorted(i) = sorted(i + 1)
             sorted(i + 1) = temp
             swap = .TRUE.
          end if
       end do
       if (.not. swap) exit
    end do

  end subroutine sort_asend

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

  subroutine create_norm_xyz(atom_counts,fil_nam,write_norm_xyz,title,ichar_title)

    character (len=100), intent(inout) :: fil_nam, write_norm_xyz
    integer, intent(in) :: atom_counts
    ! local variables
    character (len=5) :: atom_label, atom_label1
    character (len=10) :: atom_xyz_1, atom_xyz_2, atom_xyz_3
    integer ::  ierr, line_number, start_count, end_count, stat, ichars, location, location1
    integer :: geom, label_length, i, j, stats,l, structures, atom_number, ichar_title, temp_count
    character (len=11) :: search = "-----------"
    character (len=256) :: line_from_50, title
    character (len=20), dimension(:,:),allocatable :: input_coords
    real, dimension(3) :: atom_coord_1, atom_coord_2, norm_coord

! - - - - - - - - - - - Open files - - - - - - - - - - -

    label_length = len_trim(write_norm_xyz)
    atom_label = write_norm_xyz
    allocate(input_coords(4,atom_counts))
    open (status="old", unit=50, action="read",access="sequential", file="labelled.xyz", iostat=ierr)                       ! unit 50 labelled.xyz
    if ( ierr.ne.0 ) then
       write(*,'(A)') 'Unable to open labelled.xyz for creating normalised  xyz, please check file'
    else
       write(*,'(A)')' Labelled.xyz has been opened successfully.'
    end if

    stat = 0
    structures = 0
    do while (stat.eq.0)
       read(50,'(A)',iostat=stat) line_from_50
       if (stat.gt.0) then
          print*,'ERROR - something went wrong when reading labelled.xyz.'
          stop
       else if (stat.lt.0) then
          print*,'WARNING - end of labelled.xyz file reached while counting structures ensure &
               END ifthe last line of the file for labelled.xyz.'
          stop
       else
          ichars=len_trim(line_from_50)
       end if

       if (line_from_50.eq.'END') then
          print*,'All geometries cycled'
          exit
       end if

       location = index(line_from_50(1:ichars),search)

       if(location.ne.0) then
          structures = structures + 1
       end if
    end do
    rewind(50)

    open (status="replace", unit=120, action="write",access="sequential", file="NORMAL.xyz", iostat=ierr)                         ! unit 120 NORMAL.xyz
    if ( ierr.ne.0 ) then
       write(*,'(A)') 'Unable to open NORMAL.xyz for creating normalised  xyz, please check file'
    else
       write(*,'(A)')' NORMAL.xyz has been opened successfully, now creating normalised xyz file'
    end if

    open (status="replace", unit=140, action="write",access="sequential", file="MULTI.xyz", iostat=ierr)                         ! unit 140 MULTI.xyz
    if ( ierr.ne.0 ) then
       write(*,'(A)') 'Unable to open MULTI.xyz for creating normalised  xyz, please check file'
    else
       write(*,'(A)')' MULTI.xyz has been opened successfully, now creating normalised xyz file'
    end if

    open (status="replace", unit=145, action="write",access="sequential", file="NORM-MULTI.xyz", iostat=ierr)                    ! unit 144 NORM-MULTI.xyz
    if ( ierr.ne.0 ) then
       write(*,'(A)') 'Unable to open NORM-MULTI.xyz for creating normalised  xyz, please check file'
    else
       write(*,'(A)')' NORM-MULTI.xyz has been opened successfully, now creating normalised xyz file'
    end if
! - - - - - - - - - - - Initalisation - - - - - - - - - -

    line_number = 0
    geom = 0
    start_count = 0
    end_count = 0
    stats = 0
    stat = 0
    do i = 1, 3
       atom_coord_1(i) = 0.0
       atom_coord_2(i) = 0.0
       norm_coord(i) = 0.0
    end do
! - - - - - - - - - - - Read and store/write output - - -
    atom_number = atom_counts * structures
    write(120,'(I10)') atom_number
    write(120,'(A)'), title(1:ichar_title)

    do while (stat.eq.0)

       line_number = line_number + 1
       read(50,'(A)',iostat=stat) line_from_50

       if (stat.gt.0) then
          print*,'ERROR - something went wrong when reading labelled.xyz for NORMAL.xyz.'
          stop
       else if (stat.lt.0) then
          print*,'WARNING - end of labelled.xyz file reached ensure END is the last line of the file for labelled.xyz.'
          stop
       else
          ichars=len_trim(line_from_50)
       end if

       if (line_from_50(1:ichars).eq.'END') then
          print*,'All geometries cycled'
          exit
       end if

       location = index(line_from_50(1:ichars),search)

       if(location.ne.0) then
          start_count = line_number + 1
          end_count = start_count + atom_counts
          geom = geom + 1
          write(140,'(I6)') atom_counts
          write(140,'(A17,I6)') 'Snapshot number ', geom
          write(145,'(I6)') atom_counts
          write(145,'(A17,I6)') 'Snapshot number ', geom
          temp_count = start_count
          do i = 1, atom_counts
             temp_count = temp_count + 1
             read(50,'(A)',iostat=stat) line_from_50
             location1 = index(line_from_50(1:ichars),atom_label(1:label_length))
             if(location1.ne.0) then

                read(line_from_50,*) atom_label, atom_xyz_1, atom_xyz_2, atom_xyz_3

                if (index(atom_xyz_1,'E').ne.0) then
                   read(atom_xyz_1,'(e15.8)') atom_coord_1(1)
                else
                   read(atom_xyz_1,'(f15.8)') atom_coord_1(1)
                end if

                if (index(atom_xyz_2,'E').ne.0) then
                   read(atom_xyz_2,'(e15.8)') atom_coord_1(2)
                else
                   read(atom_xyz_2,'(f15.8)') atom_coord_1(2)
                end if

                if (index(atom_xyz_3,'E').ne.0) then
                   read(atom_xyz_3,'(e15.8)') atom_coord_1(3)
                else
                   read(atom_xyz_3,'(f15.8)') atom_coord_1(3)
                end if
             end if
          end do

          do i = 1, atom_counts
             backspace(50)
             temp_count = temp_count - 1
          end do
          if(start_count.ne.temp_count) then
             print*, 'Subroutine create_norm_xyz : ERROR : in locating atom coordintes to normalise against, please use &
                  MULTI.xyz as your visualisation file'
          end if
       end if

       if(line_number.ge.start_count.and.line_number.le.end_count - 1) then

          l = l + 1
          if(l.gt.atom_counts) then
             print*, 'ERROR - array assignment out of range when making NORMAL.xyz'
          end if

          read(line_from_50,*) input_coords(1,l), input_coords(2,l), input_coords(3,l), input_coords(4,l)

          read(input_coords(1,l),'(A)') atom_label1
          read(input_coords(2,l),'(f15.8)') atom_coord_2(1)
          read(input_coords(3,l),'(f15.8)') atom_coord_2(2)
          read(input_coords(4,l),'(f15.8)') atom_coord_2(3)

          do j = 1,3
             norm_coord(j) = atom_coord_2(j) - atom_coord_1(j)
          end do

          write(120,'(A3,3x,3(f15.8,3X))')atom_label1(:1), norm_coord(1), norm_coord(2), norm_coord(3)
          write(140,'(A3,3x,3(f15.8,3X))')atom_label1(:1), atom_coord_2(1), atom_coord_2(2), atom_coord_2(3)
          write(145,'(A3,3x,3(f15.8,3X))')atom_label1(:1), norm_coord(1), norm_coord(2), norm_coord(3)

          do j = 1, 3
             norm_coord(j) = 0.0
             atom_coord_2(j) = 0.0
          end do
       end if
       if(line_number.ge.end_count - 1) then
          l = 0
       end if
    end do
    close(50)
    close(120)
    close(140)
    deallocate(input_coords)
  end subroutine create_norm_xyz

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

  subroutine rdfcreate(tstep, bl, total, shell1, nx, ny, nz, solvent_atom_type, solutes_atom_type)     ! Subrountine originally created by Rob Coates Univeristy of Manchester 2015, Adapted by J. L. McDonagh
                                                                                ! R.C. calculates a rdf from an atom to surrounding water oxygen atoms
    ! Dummy variables
    integer, intent(in) :: tstep
    integer, intent(inout) :: total
    real, intent(in) :: bl
    real, intent(inout) :: shell1
    real, dimension(tstep), intent(inout) :: nx, ny, nz
    character (len=6), intent(in) :: solvent_atom_type, solutes_atom_type

    ! Local variables
    integer :: ier, i, counter, atm, stepcount, k, ncount, icharlensolv, icharlensolu
    real :: rho, dr, distance, rmax, dis, perdis
    real, dimension(:), allocatable :: x1, y1, z1
    real, dimension(:, :), allocatable :: r
    character (len=256) :: lrd
    character (len=17) :: dummy0, dummy1, dummy2, dummy3
!
! - - - - - - - - - - - - - - - - - - - - - Initialisation - - - - - - - - - - - - - - - - - - - - -
!
    ncount = 0
    k = 0
    stepcount = 0
    rmax = 0
    dr = 0
    rho = 0
    i = 0
    ier = 0
    counter = 0
    atm = 0
    shell1 = 0.0
    k = 12*bl ! divides box area by k into shells
    k = int(k)
    icharlensolv = len(trim(adjustl(solvent_atom_type)))
    icharlensolu = len(trim(adjustl(solutes_atom_type)))
!
    print*,',' , solvent_atom_type, ',' , solutes_atom_type, ','
    allocate(x1(total), y1(total), z1(total))	  ! allocates an array for storing water O coordinates at each frame

! - - - - - - - - - - - - - - - - - - - - - Radial distribution setup - - - - - - - - - - - - - - - - - -

    allocate(r(k, 4))				  ! a 2D array containing distances, probability density, 1st and second derivatives of prob density
    rmax = sqrt(3.0*((bl*0.5)**2.0))		  ! the max distance an atom can be from another in the box

    dr = rmax / k		                  ! The shell distances

    do i = 1, k					  ! Fill array with shell distances
       r(i, 1) = dr*(i-1)
       r(i, 2) = 0.0
       r(i, 3) = 0.0
       r(i, 4) = 0.0
    end do

    rho = ((total)*tstep) / bl**3	          ! Mean atom density in the box across all timesteps

! - - - - - - - - - - - Calculate the distance for each atom of type - - - - - - - - - - - - - - - - - -

    print*, '|*********** CREATING * RDF *************************************************|'
    print*, '|****************************************************************************|'
    do while (ier.eq.0)

       read(51, "(A)", iostat=ier) lrd
!       print*, ier
!       print*,lrd

       if (index(lrd, solvent_atom_type(1:icharlensolv)).gt.0) then! Store coordinate data for each water O in a frame in arrays,
          counter = counter + 1		          ! the atom label searched for may need changing depending on model used
          read(lrd,*) dummy0,x1(counter), y1(counter), z1(counter)
!          print*,counter
!          read(dummy1, "(ES17.10)") x1(counter)
!          read(dummy2, "(ES17.10)") y1(counter)
!          read(dummy3, "(ES17.10)") z1(counter)
       else if (index(lrd, solutes_atom_type(1:icharlensolu)).gt.0) then
          ncount = ncount + 1		          ! counts through the atoms, one for each frame
          read(lrd,*) dummy0, nx(ncount), ny(ncount), nz(ncount)
!          read(dummy1, "(ES17.10)") nx(ncount) ! store coordinate data for use in rdf calculation and for use in later cutting out other water molecules
!          read(dummy2, "(ES17.10)") ny(ncount)
!          read(dummy3, "(ES17.10)") nz(ncount)
       else if (index(lrd, 'Timestep').gt.0.and.stepcount.eq.0) then	! This only applies to the first line of the file
          stepcount = stepcount + 1				        ! begins counting of frames
          counter = 0
       else if (index(lrd, 'Timestep').gt.0.and.stepcount.gt.0) then	! Occurs at the end of each step

          do counter = 1, total			                        ! Calculates the distance from current target atom to each Oxygen atom in the current frame, using created function
             dis = perdist(x1(counter), y1(counter), z1(counter), nx(ncount), ny(ncount), nz(ncount), bl)

             do i = 2, k
                                                                !************|
                                                                ! Counts the |
                if (dis.ge.r(i-1,1).and.dis.lt.r(i,1)) then     !   number   |
                  r(i-1,2) = r(i-1,2) + 1                      	!     of     |
                else						                    !   Oxygen   |
                   cycle					                    !  atoms in  |
                end if						                    ! each shell |
                    					                        !************|
           end do
          end do

          counter = 0	                                                 ! reset oxygen counter for next frame

          cycle
       else
          cycle
       end if
    end do

    do i = 1, k
       r(i,2) = r(i,2) / (rho*4*pi*((r(i,1)+dr)**2)*dr)		         ! Converts the number of atoms in each shell to a probability
    end do

    do i = 2, (k-1)
       r(i,3) = r(i+1,2) - r(i-1,2) 					 ! differential
    end do

    do i = 3, k
       r(i,4) = r(i,3) - r(i-1,3) 					 ! 2nd differential
    end do

    do i = 3, k
       if (abs(r(i, 3)).lt.0.08.and.r(i-1, 3).LT.0.and.r(i, 4).GT.0) then	! this tries to find the first minimum on the graph
          shell1 = r(i, 1)						        ! May need tweaking to get correct, so check by plotting the rdf
          print*, 'First solvation shell: ', r(i, 1), ' Angstroms for atom ', solutes_atom_type  ! I've tried to make it less sensitive to small fluctuations
          exit
       end if
    end do

! - - - - - - - - - - - Write RDF data to file - - - - - - - - - - - - - - - - - - - - - - - - - -

    write(52, "(A5,12X,A5)") solutes_atom_type, solvent_atom_type
    write(52, "(A)") 'Distance    Prob. Density    Derivative    2nd Derivative'

    do i = 1, k
       write(52, "(ES17.10, 2X, ES17.10, 2X, ES17.10, 2X, ES17.10)") r(i, 1), r(i, 2), r(i, 3), r(i, 4)
    end do

    print*, '|*********** RDF * PRODUCED *************************************************|'
    print*, '|****************************************************************************|'

! - - - - - - - - - - - Clean up - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    deallocate(x1, y1, z1)

  end subroutine rdfcreate

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

  subroutine rdfcreatemicro(tstep, bl, total, shell1, nx, ny, nz, solvent_atom_type, solutes_atom_type) ! Subrountine created by Rob Coates Univeristy of Manchester 2015
    integer :: ier, i, counter, atm, total, stepcount, k, ncount, tstep
    real :: bl, rho, dr, distance, rmax, dis, shell1
    real, dimension(tstep) :: nx, ny, nz
    real, dimension(:), allocatable :: x1, y1, z1
    real, dimension(:, :), allocatable :: r
    real, parameter :: Pi = 3.1415927
    character (len=256) :: lrd
    character (len=6) :: solvent_atom_type, solutes_atom_type

! - - - - - - - - - - - - - - - - - - - - - I N I T I A L I S A T I O N - - - - - - - - - - - - - - - - - - - - -

    ncount = 0
    k = 0
    stepcount = 0
    rmax = 0
    dr = 0
    rho = 0
    i = 0
    ier = 0
    counter = 0
    atm = 0
    k = 12*int(bl) !divides box area by k into shells

! - - - - - - - - - - - - - - - - - - - - - C O U N T - A T O M - T Y P E - - - - - - - - - - - - - - - - - - - -

    allocate(x1(total), y1(total), z1(total))	!allocates an array for storing water O coordinates at each frame

! - - - - - - - - - - - - - - - - - - - - - R A D I A L - D I S - S E T U P - - - - - - - - - - - - - - - - - -

    allocate(r(k, 4))				! a 2D array containing distances, probability density, 1st and second derivatives of prob density
    rmax = sqrt(3.0*((bl*0.5)**2.0))		! the max distance an atom can be from another in the box

    dr = rmax / k					! The shell distances

    do i = 1, k					! Fill array with shell distances
       r(i, 1) = dr*(i-1)
       r(i, 2) = 0.0
       r(i, 3) = 0.0
       r(i, 4) = 0.0
    end do

    rho = ((total)*tstep) / bl**3	! Mean atom density in the box across all timesteps

! - - - - - - - - - - - C A L C U L A T E - T H E - D I S T A N C E S - F O R - E A C H - A T O M - O F - T Y P E - - - - - - -

    print*, '|*********** CREATING * RDF *************************************************|'
    print*, '|****************************************************************************|'
    do while (ier.EQ.0)

       read(51, "(A)", iostat=ier) lrd

       if (index(lrd, solvent_atom_type).GT.0) then	! Store coordinate data for each water O in a frame in arrays,
          counter = counter + 1		! the atom label searched for may need changing depending on model used

          read(lrd(7:24), "(ES17.10)") x1(counter)

          read(lrd(26:43), "(ES17.10)") y1(counter)

          read(lrd(45:62), "(ES17.10)") z1(counter)

       else if (index(lrd, solutes_atom_type).GT.0) then	!change NH3 to whichever atom is chosen as target, here Nitrogen was used
          ncount = ncount + 1		!counts through the nitrogen atoms, one for each frame

          read(lrd(7:24), "(ES17.10)") nx(ncount)		!store coordinate data for use in rdf calculation
        !and for use in later cutting out other water molecules
          read(lrd(26:43), "(ES17.10)") ny(ncount)

          read(lrd(45:62), "(ES17.10)") nz(ncount)

       else if (index(lrd, 'Timestep').GT.0.AND.stepcount.EQ.0) then	! This only applies to the first line of the file
          stepcount = stepcount + 1				! begins counting of frames
          counter = 0
       else if (index(lrd, 'Timestep').GT.0.AND.stepcount.GT.0) then	! Occurs at the end of each step

          do counter = 1, total			! Calculates the distance from current target atom
									! to each Oxygen atom in the current frame, using created function
             dis = distancemicro(x1(counter), y1(counter), z1(counter), nx(ncount), ny(ncount), nz(ncount), bl)

             do i = 1, k								       !************|
                                                               ! Counts the |
                if (dis.GE.r(i,1).AND.dis.LT.r(i+1,1)) then    !   number   |
                   r(i,2) = r(i,2) + 1                         !     of     |
                else								           !   Oxygen   |
                   cycle							           !  atoms in  |
                end if								           ! each shell |
             end do								               !************|

          end do

          counter = 0	!reset oxygen counter for next frame

          cycle
       else
          cycle
       end if
    end do

    do i = 1, k
       r(i,2) = r(i,2) / (rho*4*Pi*((r(i,1)+dr)**2)*dr)		! Converts the number of atoms in each shell to a probability
    end do

    do i = 2, (k-1)
       r(i,3) = r(i+1,2) - r(i-1,2) 					! differential
    end do

    do i = 3, k
       r(i,4) = r(i,3) - r(i-1,3) 					! 2nd differential
    end do

    do i = 3, k
       if (abs(r(i, 3)).LT.0.08.AND.r(i-1, 3).LT.0.AND.r(i, 4).GT.0) then	!this tries to find the first minimum on the graph
          shell1 = r(i, 1)						!May need tweaking to get correct, so check by plotting the rdf
          print*, 'First solvation shell: ', r(i, 1), ' Angstroms' 	!I've tried to make it less sensitive to small fluctuations
          exit
       end if
    end do

! - - - - - - - - - - - W R I T E - R D F - D A T A - T O - F I L E - - - - - - - - - - - - - - - - - - - - - - - - - -

    write(52, "(A)") "N_3			OW"
    write(52, "(A)") 'Distance    Prob. Density    Derivative    2nd Derivative'

    do i = 1, k
       write(52, "(ES17.10, 2X, ES17.10, 2X, ES17.10, 2X, ES17.10)") r(i, 1), r(i, 2), r(i, 3), r(i, 4)
    end do

    print*, '|*********** RDF * PRODUCED *************************************************|'
    print*, '|****************************************************************************|'

! - - - - - - - - - - - C L E A N - U P - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    deallocate(x1, y1, z1)

  end subroutine rdfcreatemicro

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

end module module_subroutines


text:
      SUBROUTINE rdhyd2d(numfl,headf,tailf)
      !------------------------------------------------------------------
      !Reads and initializes 2D hydroforming data (as follower load)
      !------------------------------------------------------------------
      USE loa_db, ONLY : foll_seg, add_foll
      USE c_input
      IMPLICIT NONE
      INTEGER (kind=4),INTENT(OUT):: numfl
      TYPE (foll_seg),POINTER :: headf, tailf

      !Local variables
      INTEGER (kind=4) :: chnode,i,nnode,lnofl(2)

      REAL (kind=8) :: press
      TYPE (foll_seg),POINTER :: seg

      CALL listen ('RDHYD2')
      press=getrea('PRESS ',1.d0,'!Reference pressure                ')

      CALL listen ('RDHYD2')
      IF (.NOT.exists('SURFAC')) backs = .TRUE.

      CALL listen ('RDHYD2')
      IF (exists('ELEMEN')) THEN
        nnode = 2
        NULLIFY(headf,tailf)   !Initialize follower load list

        !Loop to read surface segments
        DO
          CALL listen('RDHYD2')
          IF (exists('ENDELE')) EXIT

          ALLOCATE (seg)
          lnofl(1) = INT(param(1))
          lnofl(2) = INT(param(2))

          IF (lnofl(1)==0 .OR. lnofl(2)==0) CALL runend('RDHYD2D: END_ELEMENT_DEF expected  ')

          numfl = numfl + 1
          seg%nnode = nnode
          DO i=1,nnode
            seg%lnofl(i) = chnode(lnofl(i))
          END DO
          seg%fload = press
          CALL add_foll (seg, headf, tailf)
        END DO
        NULLIFY(seg)
      ELSE
        CALL runend('RDHYD2D: CONNECTIVITIES EXPECTED   ')
      END IF

      CALL listen ('RDHYD2')
      IF (.NOT.exists('ENDSUR')) backs = .TRUE.

      CALL listen ('RDHYD2')
      IF ( .NOT. exists('ENDHYD'))CALL runend('RDHYD2D: END_HYDROFORMING expected ')

      RETURN
      END SUBROUTINE rdhyd2d

max_stars_repo_name: jerebenitez/IFE-simpact-openfoam

max_issues_repo_name: jerebenitez/IFE-simpact-openfoam

max_forks_repo_name: jerebenitez/IFE-simpact-openfoam


text:
! { dg-do compile }
! { dg-options "-std=f2003" }
!
! PR fortran/33197
!
! Check for Fortran 2008's ATAN(Y,X) - which is equivalent
! to Fortran 77's ATAN2(Y,X).
!
real(4)    :: r4
real(8)    :: r8
complex(4) :: c4
complex(8) :: c8

r4 = atan2(r4,r4)
r8 = atan2(r8,r8)

r4 = atan(r4,r4) ! { dg-error "Too many arguments in call to 'atan'" }
r8 = atan(r8,r8) ! { dg-error "Too many arguments in call to 'atan'" }

r4 = atan2(r4,r8) ! { dg-error "same type and kind" }
r4 = atan2(r8,r4) ! { dg-error "same type and kind" }

r4 = atan2(c4,r8) ! { dg-error "must be REAL" }
r4 = atan2(c8,r4) ! { dg-error "must be REAL" }
r4 = atan2(r4,c8) ! { dg-error "same type and kind" }
r4 = atan2(r8,c4) ! { dg-error "same type and kind" }

r4 = atan2(c4,c8) ! { dg-error "must be REAL" }
r4 = atan2(c8,c4) ! { dg-error "must be REAL" }
end

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
program rungekutta
implicit none
real(kind=kind(1.0D0)) :: t,dt,tstart,tstop
real(kind=kind(1.0D0)) :: y,k1,k2,k3,k4
tstart =0.0D0 ; tstop =10.0D0 ; dt = 0.1D0
y = 1.0D0
t = tstart
write(6,'(A,f4.1,A,f12.8,A,es13.6)') 'y(',t,') = ',y,' Error = '&
           &,abs(y-(t**2+4.0d0)**2/16.0d0)
do; if ( t .ge. tstop ) exit
   k1 = f (t           , y                 )
   k2 = f (t+0.5D0 * dt, y +0.5D0 * dt * k1)
   k3 = f (t+0.5D0 * dt, y +0.5D0 * dt * k2)
   k4 = f (t+        dt, y +        dt * k3)
   y = y + dt *( k1 + 2.0D0 *( k2 + k3 ) + k4 )/6.0D0
   t = t + dt
   if(abs(real(nint(t))-t) .le. 1.0D-12) then
      write(6,'(A,f4.1,A,f12.8,A,es13.6)') 'y(',t,') = ',y,' Error = '&
           &,abs(y-(t**2+4.0d0)**2/16.0d0)
   end if
end do
contains
  function f (t,y)
    implicit none
    real(kind=kind(1.0D0)),intent(in) :: y,t
    real(kind=kind(1.0D0)) :: f
    f = t*sqrt(y)
  end function f
end program rungekutta

ext: f

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


module CommandArguments
  use KindModule
  use ConstantsModule, only: ISTDOUT, LINELENGTH, LENHUGELINE
  use VersionModule,          only: VERSION, MFVNAM, IDEVELOPMODE
  use CompilerVersion
  use SimVariablesModule,     only: simfile
  use SimModule, only: store_error, ustop, store_error_unit,                   &
                       store_error_filename
  use InputOutputModule, only: upcase
  !
  implicit none
  !
  private
  public :: GetCommandLineArguments
  !
  contains

  subroutine GetCommandLineArguments()
! ******************************************************************************
! Write information on command line arguments
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    ! -- local
    character(len=LENHUGELINE) :: line
    character(len=LINELENGTH) :: header
    character(len=LINELENGTH) :: errmsg
    character(len=LINELENGTH) :: cexe
    character(len=80) :: compiler
    character(len=20) :: cdate
    character(len=17) :: ctyp
    logical :: ltyp
    logical :: lexist
    integer(I4B) :: ipos
    integer(I4B) :: iarg
    integer(I4B) :: iterm

    integer(I4B) :: icountcmd
! ------------------------------------------------------------------------------
    !
    ! -- set mf6 executable name
    icountcmd = command_argument_count()
    call get_command_argument(0, cexe)
    cexe = adjustl(cexe)
    !
    ! -- find the program basename, not including the path (this should be
    !    mf6.exe, mf6d.exe, etc.)
    ipos = index(cexe, '/', back=.TRUE.)
    if (ipos == 0) then
      ipos = index(cexe, '\', back=.TRUE.)
    end if
    if (ipos /= 0) then
      ipos = ipos + 1
    end if
    cexe = cexe(ipos:)
    !
    ! -- write header
    call get_compile_date(cdate)
    write(header, '(a,4(1x,a),a)') &
      trim(adjustl(cexe)), '- MODFLOW',                                          &
      trim(adjustl(VERSION)), '(compiled', trim(adjustl(cdate)), ')'
    !
    ! -- set ctyp
    if (IDEVELOPMODE == 1) then
      ctyp = 'Release Candidate'
      ltyp = .TRUE.
    else
      ctyp = 'Release'
      ltyp = .FALSE.
    end if
    !
    ! -- Read remaining arguments
    iarg = 0
    iterm = 0
    do iarg = 1, icountcmd
      call get_command_argument(iarg, line)
      call upcase(line)
      iterm = 1
      select case(trim(adjustl(line)))
        case('-H', '-?', '--HELP')
          call write_usage(trim(adjustl(header)), trim(adjustl(cexe)))
        case('-V', '--VERSION')
          write(ISTDOUT,'(2a,2(1x,a))') &
            trim(adjustl(cexe)), ':', trim(adjustl(VERSION)), ctyp
        case('-DEV', '--DEVELOP')
          write(ISTDOUT,'(2a,g0)') &
            trim(adjustl(cexe)), ': develop version ', ltyp
        case('-C', '--COMPILER')
          call get_compiler(compiler)
          write(ISTDOUT,'(2a,1x,a)') &
            trim(adjustl(cexe)), ':', trim(adjustl(compiler))
        case default
          call write_usage(trim(adjustl(header)), trim(adjustl(cexe)))
          write(errmsg, '(2a,1x,a)') &
            trim(adjustl(cexe)), ': illegal option -', trim(adjustl(line))
          call store_error(errmsg)
      end select
    end do
    !
    ! -- no command line arguments - check if mfsim.nam exists
    if  (icountcmd == 0) then
      inquire(file=simfile, exist=lexist)
      if (.NOT. lexist) then
        iterm = 1
        write(errmsg, '(2a,2(1x,a))') &
           trim(adjustl(cexe)), ':', simfile, 'is not present in working directory.'
        call store_error(errmsg)
      end if
    end if
    !
    ! -- command line arguments present or mfsim.nam file does not exist
    if (iterm > 0) then
      call USTOP()
    end if
    !
    ! -- return
    return
  end subroutine GetCommandLineArguments

  subroutine write_usage(header, cexe)
    ! -- dummy
    character(len=*), intent(in) :: header
    character(len=*), intent(in) :: cexe
    ! -- local
    character(len=*), parameter :: OPTIONSFMT =                                    &
      "(/,                                                                      &
      &'Options   GNU long option   Meaning ',/,                                &
      &' -h, -?   --help            Show this message',/,                       &
      &' -v       --version         Display program version information.',/,    &
      &' -dev     --develop         Display program develop option mode.',/,    &
      &' -c       --compiler        Display compiler information.',/,           &
      &'                                                                    ',/,&
      &'Bug reporting and contributions are welcome from the community. ',/,    &
      &'Questions can be asked on the issues page[1]. Before creating a new',/, &
      &'issue, please take a moment to search and make sure a similar issue',/, &
      &'does not already exist. If one does exist, you can comment (most',/,    &
      &'simply even with just :+1:) to show your support for that issue.',/,    &
      &'                                                                    ',/,&
      &'[1] https://github.com/MODFLOW-USGS/modflow6/issues',/)"
! ------------------------------------------------------------------------------
    write(ISTDOUT,'(a,/,a,1x,a,15x,a,2(1x,a),2a,/,a,1x,a,1x,a,5x,a)') &
      trim(adjustl(header)),                                              &
      'usage:', cexe, 'run MODFLOW', trim(adjustl(MFVNAM)),               &
      'using "', trim(adjustl(simfile)), '"',                             &
      '   or:', cexe, '[options]',                                        &
      'retrieve program information'
    write(ISTDOUT, OPTIONSFMT)
    !
    ! -- return
    return
  end subroutine write_usage

end module CommandArguments


text:
# Suggested simulator command:
#    VCS:
#       vcs -full64 -sverilog -ntb_opts uvm-1.2 -R -f flist.f +UVM_TESTNAME=minitb_dummy_test
#    Irun/Xcelium:
#       irun -64bit -uvm -uvmhome CDNS-1.2 $UVM_HOME/src/uvm_macros.svh -f flist.f +UVM_TESTNAME=minitb_dummy_test
#

-top minitb_tb_top

+incdir+dv/common/data/seq_lib
+incdir+dv/common/data/agent
+incdir+dv/minitb/common/vseqr
+incdir+dv/minitb/common/env
+incdir+dv/minitb/seq_lib
+incdir+dv/minitb/test_lib

rtl/minitb/minitb.sv
rtl/minitb/minitb_top.sv
dv/common/data/data_if.sv
dv/common/data/seq_lib/data_seq_lib.sv
dv/common/data/agent/data_agent.sv
dv/minitb/common/vseqr/minitb_vseqr.sv
dv/minitb/common/env/minitb_env.sv
dv/minitb/seq_lib/minitb_seq_lib.sv
dv/minitb/test_lib/minitb_test_lib.sv
dv/minitb/minitb_tb_top.sv

max_stars_repo_name: verilogtools/minitb

max_issues_repo_name: verilogtools/minitb

max_forks_repo_name: verilogtools/minitb


text:
      subroutine preadstr(word,latt,twiss,mult,istr,nstr,lfno)
      use tfstk
      use ffs
      use tffitcode
      use tfcsi,only:cssetp
      logical lod,lex,stat,name
      character*(*) word
      character*80 outc,ename,en
      character*256 tfconvstr
C*DEC
      character*256 word2
C*DEC
      integer*8 latt(nlat)
      dimension twiss(nlat,-ndim:ndim,ntwissfun),mult(*),
     &          istr(nstra,4)
      external pack
      include 'inc/common.inc'
      include 'inc/coroper.inc'
      data in/50/
c
c     call getwdlK(word)
c     call getwdl(word)
      itype=itfpeeko(ia,x,next)
      lfni1=x+.5d0
      if(itype.eq.1) then
        call cssetp(next)
        write(word,'(''ftn'',i2.2)') lfni1
      elseif(itype.eq.101) then
        call cssetp(next)
        word=tfconvstr(101,ia,x,nc,'*')
        if(word.eq.' ') then
          call permes('?Missing filename for MONREAD','.',' ',lfno)
          call getwdl(word)
          return
        endif
        call texpfn(word)
      endif
      inquire(unit=in,opened=lod,iostat=ios)
      if(lod) close(in)
C*DEC
      word2 = word
      inquire(iostat=ios,file=word2,exist=lex)
C*DEC
C*HP
C     inquire(iostat=ios,file=word,exist=lex)
C*HP
      if(.not.lex) then

        call permes('?File',word,'not found.',lfno)
        call getwdl(word)
        return
      endif
c
      open (in,file=word,iostat=ios,status='OLD')
      if(ios .ne. 0) then
        call permes('?Cannot open',word,'.',lfno)
        call getwdl(word)
        return
      endif
c
c     .... reset current specification of Coorectors
      do 10 i=1,nstra
        istr(istr(i,2),3)=1
 10   continue
      stat=.true.
      name=.true.
      jm=0
 1    continue
      call preabuf(outc,in,stat)
      if(stat) goto 99
      if(name) then
        ename=outc
      else
        read(outc,*,err=98) va1
      endif
      name=.not.name
      if(.not.name) goto 1
      ls=jm
      do 20 i=ls+1,nstra
        call elnameK(latt,istr(istr(i,2),1),mult,en)
        if(ename.eq.en) then
          jm=i
          istr(istr(i,2),3)=0
          goto 22
        endif
 20   continue
      do 21 i=ls,1,-1
        call elnameK(latt,istr(istr(i,2),1),mult,en)
        if(ename.eq.en) then
          jm=i
          istr(istr(i,2),3)=0
          goto 22
        endif
 21   continue
      call permes('??',ename,' is not a BEND.',lfno)
      goto 1
c
 22   continue
      rlist(latt(istr(istr(jm,2),1))+11)=va1
      goto 1
 98   call permes('??',outc,' is not a number.',lfno)
 99   call pack(istr(1,2),istr(1,3),nstr,nstra)
      write(lfno,'(2(A,I4))')
     1  '  Correction dipole available :',nstr,' in ',nstra
      nstract=nstr
c ....write data in the OPERation buffer.
      if(nstr.ne.0) then
        nstope=nstr
        if(istope.ne.0) then
          call tfree(int8(istope))
        endif
        istope=italoc((nstr+1)/2)
        do 120 i=1,nstr
          ilist(mod(i-1,2)+1,istope+(i-1)/2)=istr(i,2)
 120    continue
      endif
      call getwdl(word)
      close(in)
      return
      end

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD


text:
Bernards Tire, Brake & Alignment began as a familyrun operation in 1973. The shop was sold in 2006 when the owner retired.

Bernards provides general auto repair, but specializes in repairs to your vehicles drive system. Services include repair and replacement of drive system, wheels, wheel alignment, retreads, and performance rims. Much of Bernards business is by referral, and oftentimes comes by way of Redfield Machine. City of Davis trucks have also been seen serviced there. Bernards offers reasonable prices/estimates and good customer service.

The shop carries Cooper tires almost exclusively, and mounting is not included in the quoted tire price; however, Bernards includes a lifetime of free tire rotations every 5,000 miles. Cooper tires are rated highly in http://www.jdpower.com/corporate/news/releases/pressrelease.asp?ID2002050 consumer satisfaction by JD Powers.

To find more automobile repair locations within Davis, please visit the Automobile Repair page.

20051130 22:13:46 nbsp Ive found the owner to be honest and competent.  Some of the staff arent so impressive, but the owner is the one who really matters. Users/GrahamFreeman

20051220 20:37:43 nbsp Ive bought tires here three or four times over the years  always got a good price, good work done. and all in a timely manner. Highly recommended. Users/PaulThober

20060103 21:05:02 nbsp My whole family has gone to Bernards for countless years and weve always left happy.  Customer Service A+.  Seeing this reminds me I should probably get my tires rotated soon and the brakes checked!  Though I havent done it, they do offer to give you a ride to work/home after dropping your car off.  I highly recommend that everyone go to Bernards.   Users/BrettCoder

20060501 18:40:27 nbsp They try to explain everything and they are really honest. Users/AmirArdeshir

20060524 08:41:24 nbsp I just had some new tires put on my car, and am very happy with the service I received at Bernards.  They have very competitive prices (even when compared to WalMart or Costco).  They seem to sell mostly Cooper tires, but they sell Cooper offbrands as well.  When you buy tires at a place like this, you know your not going to get some cheap crap that nobody will be accountable for.  When you need tires, I highly recommend getting a quote from Bernards  . . . youll probably be surprised.   Users/AdamSchneider

20060929 21:32:23 nbsp Took my 1999 SUV in because the ABS light was on (I found them here on DavisWiki). They fixed the problem and had the car back to me the same day. They took the time to explain what was going wrong and what they did to fix the problem. I was very pleased with the service. Users/DanielBoxwell

20061101 12:46:38 nbsp They got me in for an alignmetn first thing this morning on very short notice.  The technician took the time to get the steering wheel perfectly straight, even though that meant putting the car back on the rack after his test drive.  They also had no obvious problems getting my lowerthanstock car on the rack.  Very reasonably priced, as well!   Users/RussBowlus

20061206 03:52:18 nbsp Got a 4 wheel alignment after finally chewing through a pair of tires.  Afterwards the car acted oddly (severe diving to the left).  Brought it back and they threw it up on a lift and found the right ball joint had come loose.  Fixed it, realigned and refused to take any payment for the second round.  These guys are worth doing business with. Users/TarZxf

20070228 10:19:48 nbsp These guys go above and beyond.  They came to pick up my car for me while I was at school, and saved me a lot of trouble.  Higly recommended! Users/BeckiBrunelli

20070827 09:38:07 nbsp I had my tires taken care of by Bernards for years, and I was always pleased with the product and service. The owner, Ed, took over the business from his dad years ago. I learned last week that Ed retired and sold the place last year. Some of the same guys work there, but I can no longer endorse the shop because its under new ownership. Users/robinlaughlin

20071107 14:07:12 nbsp Called today to check prices on tire rotation and balancing.  Bernards offers tire rotation and balancing of the front tires for $28.50.  Additional tires are balanced for $12.95/tire. Users/MattVanGundy

20080303 22:01:31 nbsp They seem fine. They realigned my tires and fixed the rattling. Users/JohnWong

20080820 14:12:29 nbsp Good people to do business with.  Clear prices, service when promised and often sooner.  When my 6 month old tires from them started going flat frequently, they cleaned up and sealed the rims for free, no questions (it was my old rims, not their new tires at fault).   Users/NotTires

20081115 14:59:31 nbsp I put in new tie rod ends last night and threw my alignment way out of whack (as expected). Dropped by Bernards with no appointment and they got my truck aligned in about an hour and a half, and it tracks perfectly straight. I also asked if they could grease my new tie rod ends and they did that at no extra charge. Very nice guys there sat and BSd about tires with one of the employees while another was greasing the tie rods. Very friendly. Will do business with them again. Users/ChrisC

20081115 19:51:46 nbsp Very friendly customer service and quick turnaround. I bought a set of new tires and got an alignment. They gave me a quote for 3 different tires and went over the merits of each set. I was very pleased with their service. Users/evergeneva

20090315 11:14:55 nbsp I like this place. Theyre understanding and willing to go the distance to make sure your vehicle is fixed with as little hassle possible. Ill go back again.  Users/Jeserah

20090407 17:12:08 nbsp GREAT experience! I was referred by some friends and I was given a free brake inspection (promo.) When I went to pick up my truck, I learned that I still had 70% brake life left and I dont need to repair anything, whereas a $15 brake inspection at Sears got me a quote saying that I had less than 25% brake life left. Just goes to show the honesty of this shop and that theyre not trying to squeeze money out of you by saying you need to replace things that dont need to be, unlike OTHER places. When my brakes DO wear out, I know where IM going for honest service! Thanks Bernards!!! Users/Omicron

20090601 20:15:24 nbsp While the manager Yusef seems nice, I have to question the integrity of the mechanics. A routine wheel bearing job suspiciously lead to a failed transmission moments after receiving the car back from them. Also, after dealing with the transmission fiasco, on the same wheel the bearing was replaced, the brakes are now squeaking. I had asked them to check it out for me since they worked on the components last, but the manager was quick to blame the brakes, even though they had no problem before. After bringing the car to my normal mechanic, it turns out the caliper pins in the brakes were not greased by Bernards mechanic and was sticking. (The brakes have to come off to do a wheel bearing). Would I go back here again? No. Would I purchase tires? Maybe. If you have a european import, I suggest staying away from this place. They are good for tires and alignments, but I would stay away for any mechanical jobs. Users/lawrenceucd

20092508 23:15:24 nbsp STAY AWAY from this place. Had to do a regular check (no damage nor break) and asked for the approx. price. They told me between 150220$ which seamed reasonable (checked some adds before). Let the car there and told them to call me, in case they found smth. unexpected. Called me to tell me that all was right (great car!) and they dont even had to change the plugs as they were still good (great service, I thought). They also told me the price, but couldnt fully understand on phone (foreigner). Sounded as in the expected range. However, when I came in to pick up the car the costs where 5 TIMES higher as previously told! Wouldnt pay and resisted to the told price ... however, they said they had told me the price on phone so had no chance but to pay. They gave me a 10% discount, but I still felt really ripped OFF. The mechanic said he had been working the whole day on the car (to change oil of engine, differencials, and for general service ... OF COURSE!!!). For those of you who still want to go there > the rip seams to work as follows: Tell the car owner that it will cost little, show the owner that you are willing to save his money by not doing all the required changes, mumble the price during phone call, say you have worked the whole day on the car because with these new cars everything takes much longer, add high oil prices and say you dont get better prices yourself, and HAVE A GREAT NIGHT OUT after the ripped owner walked from yard. STAY AWAY from this place! Users/tfb676

20091111 11:03:21 nbsp Seems like a nice place up front until you find out that they charge you almost $80 just to look at your car.  I know thats a huge red flag right there, but I decided to have them do the labor anyway based on all the good reviews theyve received on here.  Also a poor choice.  The labor costs at Bernards for a simple break job are nearly $90/hr.  I payed $437.25 in labor on a job that took on 4 hours to do.  Seemed to have done a good job, but not worth the cost.  Take your business to the honest guy wholl take a look and give you an estimate for free. Users/LoganTurner

20091228 18:27:54 nbsp I am really impressed with Bernards! I called around to compare prices and try to find the best deal.  Not only were they the most competitive, they saved me from buying tires for $500+ that are only good for ~ 25,000 miles.  The technician at Goodyear told me what a great deal the tires were and I was going to buy them.  When I talked to Bernards, he spent a good amount of time going through their inventory trying to find an ideal tire for me.  When I asked the difference between the tire he recommend and the good deal tire for Goodyear, I was shocked when he told me the life span was half the other tire!  I feel like the Goodyear worker was not being honest.  I am VERY HAPPY I called Bernards.    Users/CharlieMurphy

20091228 19:53:54 nbsp How many Bernards are there? Users/CovertProfessor

20100715 13:29:39 nbsp Good prices and quick service but dont leave any change in your car or it will be gone, EVERY TIME!! Users/rdgalh@att.net

20110223 10:05:12 nbsp They installed tires for me. The whole process was fast and professional. They did not charge much  I thought price was fair. Also, I asked to fix headlight cords that were dangling, they did it for free. I will be coming back for sure. Thank you Bernards!  Users/Adelina

20110809 19:10:35 nbsp My brakes started feeling funny on my truck and after researching it on the internet, I found that it was likely to be my brake booster. Rather then change it out and find out it was some other part, I looked for a locally owned shop in Davis to verify if it was indeed the brake booster or not.  Bernards seemed to have positive reviews so I decided to go there.  They gave a free estimate and told me that it was a check valve on the brake booster, not the actual brake booster. The part was $40 dollars and they charged me $47.50 to install it. My brakes started acting up again about 2 days later. I took it back to Bernards and they then told me it was indeed my brake booster that was bad and it would cost me $550 dollars to fix it. I decided to wait and went home and watched a youtube video on how to change a brake booster. The extent of the labor was removing 8 bolts, switching out the brake booster, and then reinstalling the 8 bolts. I went to a used auto parts place in Sacramento and got a brake booster for $75 dollars. When I went to look under my hood, I saw the check valve Bernards had previously installed had 2 metal clips on it. They had charged me $47.50 in labor to remove two clips and put the check valve in! This is something you can do using your fingers. I ended up installing the brake booster in 2 hours using my cheap tools I bought from Kragens auto parts. Since I changed the brake booster in 2 hours, I dont understand how Bernards can justify charging $550 dollars to change it using their skilled mechanics and nice tools. They can probably change it in 45 minutes if even that.  I regret that they got my money to change the check valve that did not actually need to be replaced , however, the experience did give me a real good impression of the things going on there. I prefer to support independent businesses owned by local people but Bernards is one place that will never get my recommendation or my money again. I suggest that you look elsewhere for a local shop to fix your car.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


!! Thoroughly test send, i.e. foo[N] = bar in all variants
!!
!! Do simple tests for send(). These test comprise
!!
!! FOO [N] = BAR
!!
!! where
!!
!!  FOO                BAR               images
!! scalar            scalar             N == me
!!  int(k e [1,4])    int(k e [1,4])
!!  real(k e [4,8])   real(k e [4,8])
!!  int(k e [1,4])    real(k e [4,8])
!!  real(k e [4,8])   int(k e [1,4])
!!
!! array(1:5)        scalar
!!  int(k e [1,4])    int(k e [1,4])
!!  real(k e [4,8])   real(k e [4,8])
!!  int(k e [1,4])    real(k e [4,8])
!!  real(k e [4,8])   int(k e [1,4])
!!
!! array(1:5)        array(1:5)
!!  int(k e [1,4])    int(k e [1,4])
!!  real(k e [4,8])   real(k e [4,8])
!!  int(k e [1,4])    real(k e [4,8])
!!  real(k e [4,8])   int(k e [1,4])
!!
!! array(1:3)       array(::2)
!!  int(k e [1,4])    int(k e [1,4])
!!  real(k e [4,8])   real(k e [4,8])
!!  int(k e [1,4])    real(k e [4,8])
!!  real(k e [4,8])   int(k e [1,4])
!!
!! array(4:5)       array(2::2)
!!  int(k e [1,4])    int(k e [1,4])
!!  real(k e [4,8])   real(k e [4,8])
!!  int(k e [1,4])    real(k e [4,8])
!!  real(k e [4,8])   int(k e [1,4])
!!
!! array(1:3)      array(3:1:-1)
!!  int(k e [1,4])    int(k e [1,4])
!!  real(k e [4,8])   real(k e [4,8])
!!  int(k e [1,4])    real(k e [4,8])
!!  real(k e [4,8])   int(k e [1,4])
!!
!! all of the above but for              N != me
!!
!! And may be some other, I've forgotten.
!!
!! Author: Andre Vehreschild, 2017

program send_convert_nums

  implicit none

  real(kind=4), parameter :: tolerance4 = 1.0e-4
  real(kind=8), parameter :: tolerance4to8 = 1.0E-4
  real(kind=8), parameter :: tolerance8 = 1.0E-6

  integer(kind=1), allocatable, codimension[:] :: co_int_scal_k1
  integer(kind=1), allocatable                 :: int_scal_k1
  integer(kind=4), allocatable, codimension[:] :: co_int_scal_k4
  integer(kind=4), allocatable                 :: int_scal_k4
  real(kind=4)   , allocatable, codimension[:] :: co_real_scal_k4
  real(kind=4)   , allocatable                 :: real_scal_k4
  real(kind=8)   , allocatable, codimension[:] :: co_real_scal_k8
  real(kind=8)   , allocatable                 :: real_scal_k8

  integer(kind=1), allocatable, dimension(:), codimension[:] :: co_int_k1
  integer(kind=1), allocatable, dimension(:)                 :: int_k1
  integer(kind=4), allocatable, dimension(:), codimension[:] :: co_int_k4
  integer(kind=4), allocatable, dimension(:)                 :: int_k4
  real(kind=4)   , allocatable, dimension(:), codimension[:] :: co_real_k4
  real(kind=4)   , allocatable, dimension(:)                 :: real_k4
  real(kind=8)   , allocatable, dimension(:), codimension[:] :: co_real_k8
  real(kind=8)   , allocatable, dimension(:)                 :: real_k8

  logical :: error_printed=.false.

  associate(me => this_image(), np => num_images())
    if (np < 2) error stop 'Can not run with less than 2 images.'

    allocate(int_scal_k1, SOURCE=INT(42, 1))
    allocate(int_scal_k4, SOURCE=42)
    allocate(co_int_scal_k1[*]) ! allocate syncs here
    allocate(co_int_scal_k4[*]) ! allocate syncs here
    allocate(int_k1(1:5), SOURCE=INT([ 5, 4, 3, 2, 1], 1))
    allocate(int_k4(1:5), SOURCE=[ 5, 4, 3, 2, 1])
    allocate(co_int_k1(5)[*]) ! allocate syncs here
    allocate(co_int_k4(5)[*]) ! allocate syncs here

    allocate(real_scal_k4, SOURCE=37.042)
    allocate(real_scal_k8, SOURCE=REAL(37.042, 8))
    allocate(co_real_scal_k4[*]) ! allocate syncs here
    allocate(co_real_scal_k8[*]) ! allocate syncs here
    allocate(real_k4(1:5), SOURCE=[ 5.1, 4.2, 3.3, 2.4, 1.5])
    allocate(real_k8(1:5), SOURCE=REAL([ 5.1, 4.2, 3.3, 2.4, 1.5], 8))
    allocate(co_real_k4(5)[*]) ! allocate syncs here
    allocate(co_real_k8(5)[*]) ! allocate syncs here
    ! First check send/copy to self
    if (me == 1) then
      co_int_scal_k1[1] = int_scal_k1
      print *, co_int_scal_k1
      if (co_int_scal_k1 /= int_scal_k1) call print_and_register( 'send scalar int kind=1 to kind=1 self failed.')

      co_int_scal_k4[1] = int_scal_k4
      print *, co_int_scal_k4
      if (co_int_scal_k4 /= int_scal_k4) call print_and_register( 'send scalar int kind=4 to kind=4 self failed.')

      co_int_scal_k4[1] = int_scal_k1
      print *, co_int_scal_k4
      if (co_int_scal_k4 /= int_scal_k4) call print_and_register( 'send scalar int kind=1 to kind=4 self failed.')

      co_int_scal_k1[1] = int_scal_k4
      print *, co_int_scal_k1
      if (co_int_scal_k1 /= int_scal_k1) call print_and_register( 'send scalar int kind=4 to kind=1 self failed.')

      co_int_k1(:)[1] = int_k1
      print *, co_int_k1
      if (any(co_int_k1 /= int_k1)) call print_and_register( 'send int kind=1 to kind=1 self failed.')

      co_int_k4(:)[1] = int_k4
      print *, co_int_k4
      if (any(co_int_k4 /= int_k4)) call print_and_register( 'send int kind=4 to kind=4 self failed.')

      co_int_k4(:)[1] = int_k1
      print *, co_int_k4
      if (any(co_int_k4 /= int_k4)) call print_and_register( 'send int kind=1 to kind=4 self failed.')

      co_int_k1(:)[1] = int_k4
      print *, co_int_k1
      if (any(co_int_k1 /= int_k1)) call print_and_register( 'send int kind=4 to kind=1 self failed.')
    else if (me == 2) then ! Do the real copy to self checks on image 2
      co_real_scal_k4[2] = real_scal_k4
      print *, co_real_scal_k4
      if (abs(co_real_scal_k4 - real_scal_k4) > tolerance4) &
        call print_and_register( 'send scalar real kind=4 to kind=4 self failed.')

      co_real_scal_k8[2] = real_scal_k8
      print *, co_real_scal_k8
      if (abs(co_real_scal_k8 - real_scal_k8) > tolerance8) &
        call print_and_register( 'send scalar real kind=8 to kind=8 self failed.')

      co_real_scal_k8[2] = real_scal_k4
      print *, co_real_scal_k8
      if (abs(co_real_scal_k8 - real_scal_k8) > tolerance4to8) &
        call print_and_register( 'send scalar real kind=4 to kind=8 self failed.')

      co_real_scal_k4[2] = real_scal_k8
      print *, co_real_scal_k4
      if (abs(co_real_scal_k4 - real_scal_k4) > tolerance4) &
        call print_and_register( 'send scalar real kind=8 to kind=4 self failed.')

      co_real_k4(:)[2] = real_k4
      print *, co_real_k4
      if (any(abs(co_real_k4 - real_k4) > tolerance4)) &
        call print_and_register( 'send real kind=4 to kind=4 self failed.')

      co_real_k8(:)[2] = real_k8
      print *, co_real_k8
      if (any(abs(co_real_k8 - real_k8) > tolerance8)) call print_and_register( 'send real kind=8 to kind=8 self failed.')

      co_real_k8(:)[2] = real_k4
      print *, co_real_k8
      if (any(abs(co_real_k8 - real_k8) > tolerance4to8)) call print_and_register( 'send real kind=4 to kind=8 self failed.')

      co_real_k4(:)[2] = real_k8
      print *, co_real_k4
      if (any(abs(co_real_k4 - real_k4) > tolerance4)) call print_and_register( 'send real kind=8 to kind=4 self failed.')
    end if

    sync all
    if (me == 1) then
      co_int_scal_k1[2] = int_scal_k1

      co_int_scal_k4[2] = int_scal_k4

      co_int_k1(:)[2] = int_k1

      co_int_k4(:)[2] = int_k4

      co_real_scal_k4[2] = real_scal_k4

      co_real_scal_k8[2] = real_scal_k8

      co_real_k4(:)[2] = real_k4

      co_real_k8(:)[2] = real_k8
    end if

    sync all
    if (me == 2) then
      print *, co_int_scal_k1
      if (co_int_scal_k1 /= int_scal_k1) call print_and_register( 'send scalar int kind=1 to kind=1 to image 2 failed.')

      print *, co_int_scal_k4
      if (co_int_scal_k4 /= int_scal_k4) call print_and_register( 'send scalar int kind=4 to kind=4 to image 2 failed.')

      print *, co_int_k1
      if (any(co_int_k1 /= int_k1)) call print_and_register( 'send int kind=1 to kind=1 to image 2 failed.')

      print *, co_int_k4
      if (any(co_int_k4 /= int_k4)) call print_and_register( 'send int kind=4 to kind=4 to image 2 failed.')

      print *, co_real_scal_k4
      if (abs(co_real_scal_k4 - real_scal_k4) > tolerance4) &
        call print_and_register( 'send scalar real kind=4 to kind=4 to image 2 failed.')

      print *, co_real_scal_k8
      if (abs(co_real_scal_k8 - real_scal_k8) > tolerance8) &
        call print_and_register( 'send scalar real kind=8 to kind=8 to image 2 failed.')

      print *, co_real_k4
      if (any(abs(co_real_k4 - real_k4) > tolerance4)) call print_and_register( 'send real kind=4 to kind=4 to image 2 failed.')

      print *, co_real_k8
      if (any(abs(co_real_k8 - real_k8) > tolerance8)) call print_and_register( 'send real kind=8 to kind=8 to image 2 failed.')
    end if

    sync all
    if (me == 1) then
      co_int_scal_k4[2] = int_scal_k1

      co_int_scal_k1[2] = int_scal_k4

      co_int_k4(:)[2] = int_k1

      co_int_k1(:)[2] = int_k4

      co_real_scal_k8[2] = real_scal_k4

      co_real_scal_k4[2] = real_scal_k8

      co_real_k8(:)[2] = real_k4

      co_real_k4(:)[2] = real_k8
    end if

    sync all
    if (me == 2) then
      print *, co_int_scal_k4
      if (co_int_scal_k4 /= int_scal_k4) call print_and_register( 'send scalar int kind=1 to kind=4 to image 2 failed.')

      print *, co_int_scal_k1
      if (co_int_scal_k1 /= int_scal_k1) call print_and_register( 'send scalar int kind=4 to kind=1 to image 2 failed.')

      print *, co_int_k4
      if (any(co_int_k4 /= int_k4)) call print_and_register( 'send int kind=1 to kind=4 to image 2 failed.')

      print *, co_int_k1
      if (any(co_int_k1 /= int_k1)) call print_and_register( 'send int kind=4 to kind=1 to image 2 failed.')

      print *, co_real_scal_k8
      if (abs(co_real_scal_k8 - real_scal_k8) > tolerance4to8) &
        call print_and_register( 'send scalar real kind=4 to kind=8 to image 2 failed.')

      print *, co_real_scal_k4
      if (abs(co_real_scal_k4 - real_scal_k4) > tolerance4) &
        call print_and_register( 'send scalar real kind=8 to kind=4 to image 2 failed.')

      print *, co_real_k8
      if (any(abs(co_real_k8 - real_k8) > tolerance4to8)) &
        call print_and_register( 'send real kind=4 to kind=8 to image 2 failed.')

      print *, co_real_k4
      if (any(abs(co_real_k4 - real_k4) > tolerance4)) call print_and_register( 'send real kind=8 to kind=4 to image 2 failed.')
    end if

    ! Scalar to array replication
    sync all
    if (me == 1) then
      co_int_k4(:)[2] = int_scal_k4

      co_int_k1(:)[2] = int_scal_k1

      co_real_k8(:)[2] = real_scal_k8

      co_real_k4(:)[2] = real_scal_k4
    end if

    sync all
    if (me == 2) then
      print *, co_int_k4
      if (any(co_int_k4 /= int_scal_k4)) call print_and_register( 'send int scal kind=4 to array kind=4 to image 2 failed.')

      print *, co_int_k1
      if (any(co_int_k1 /= int_scal_k1)) call print_and_register( 'send int scal kind=1 to array kind=1 to image 2 failed.')

      print *, co_real_k8
      if (any(abs(co_real_k8 - real_scal_k8) > tolerance8)) &
        call print_and_register( 'send real kind=8 to array kind=8 to image 2 failed.')

      print *, co_real_k4
      if (any(abs(co_real_k4 - real_scal_k4) > tolerance4)) &
        call print_and_register( 'send real kind=4 to array kind=4 to image 2 failed.')
    end if

    ! and with kind conversion
    sync all
    if (me == 1) then
      co_int_k4(:)[2] = int_scal_k1

      co_int_k1(:)[2] = int_scal_k4

      co_real_k8(:)[2] = real_scal_k4

      co_real_k4(:)[2] = real_scal_k8
    end if

    sync all
    if (me == 2) then
      print *, co_int_k4
      if (any(co_int_k4 /= int_scal_k4)) call print_and_register( 'send int scal kind=1 to array kind=4 to image 2 failed.')

      print *, co_int_k1
      if (any(co_int_k1 /= int_scal_k1)) call print_and_register( 'send int scal kind=4 to array kind=1 to image 2 failed.')

      print *, co_real_k8
      if (any(abs(co_real_k8 - real_scal_k8) > tolerance8)) &
        call print_and_register( 'send real kind=4 to array kind=8 to image 2 failed.')

      print *, co_real_k4
      if (any(abs(co_real_k4 - real_scal_k4) > tolerance4)) &
        call print_and_register( 'send real kind=8 to array kind=4 to image 2 failed.')
    end if
    ! and with type conversion
    sync all
    if (me == 1) then
      co_int_k4(:)[2] = real_scal_k4

      co_int_k1(:)[2] = real_scal_k4

      co_real_k8(:)[2] = int_scal_k4

      co_real_k4(:)[2] = int_scal_k4
    end if

    sync all
    if (me == 2) then
      print *, co_int_k4
      if (any(co_int_k4 /= INT(real_scal_k4, 4))) &
        call print_and_register( 'send real scal kind=4 to int array kind=4 to image 2 failed.')

      print *, co_int_k1
      if (any(co_int_k1 /= INT(real_scal_k4, 1))) &
        call print_and_register( 'send real scal kind=1 to int array kind=1 to image 2 failed.')

      print *, co_real_k8
      if (any(abs(co_real_k8 - int_scal_k4) > tolerance4to8)) &
        call print_and_register( 'send int kind=4 to real array kind=8 to image 2 failed.')

      print *, co_real_k4
      if (any(abs(co_real_k4 - int_scal_k4) > tolerance4)) &
        call print_and_register( 'send int kind=4 to real array kind=4 to image 2 failed.')
    end if

    ! Now with strides
    ! First check send/copy to self
    if (me == 1) then
      co_int_k1 = -1
      co_int_k1(::2)[1] = int_k1(1:3)
      print *, co_int_k1
      if (any(co_int_k1 /= [int_k1(1), INT(-1, 1), int_k1(2), INT(-1, 1), int_k1(3)])) &
        & call print_and_register( 'send strided int kind=1 to kind=1 self failed.')

      co_int_k4 = -1
      co_int_k4(::2)[1] = int_k4
      print *, co_int_k4
      if (any(co_int_k4 /= [int_k4(1), -1, int_k4(2), -1, int_k4(3)])) &
        call print_and_register( 'send strided int kind=4 to kind=4 self failed.')

      co_int_k4 = -2
      co_int_k4(2::2)[1] = int_k1(4:5)
      print *, co_int_k4
      if (any(co_int_k4 /= [-2, int_k4(4), -2, int_k4(5), -2])) &
        call print_and_register( 'send strided int kind=1 to kind=4 self failed.')

      co_int_k1 = -2
      co_int_k1(2::2)[1] = int_k4(4:5)
      print *, co_int_k1
      if (any(co_int_k1 /= [INT(-2, 1), int_k1(4), INT(-2, 1), int_k1(5), INT(-2, 1)])) &
        & call print_and_register( 'send strided int kind=4 to kind=1 self failed.')

      co_int_k1(1:5) = int_k1(5:1:-1)
      co_int_k1(::2)[1] = co_int_k1(3:1:-1)
      print *, co_int_k1
      ! Note, indezes two times reversed!
      if (any(co_int_k1 /= [int_k1(3), int_k1(4), int_k1(4), int_k1(2), int_k1(5)])) &
        & call print_and_register( 'send strided with temp int kind=1 to kind=1 self failed.')

      co_int_k4(1:5) = int_k4(5:1:-1)
      co_int_k4(::2)[1] = co_int_k4(3:1:-1)
      print *, co_int_k4
      if (any(co_int_k4 /= [int_k4(3), int_k4(4), int_k4(4), int_k4(2), int_k4(5)])) &
       & call print_and_register( 'send strided with temp int kind=4 to kind=4 self failed.')
    else if (me == 2) then ! Do the real copy to self checks on image 2
      co_real_k4 = -1.0
      co_real_k4(::2)[2] = real_k4(1:3)
      print *, co_real_k4
      if (any(abs(co_real_k4 - [real_k4(1), -1.0, real_k4(2), -1.0, real_k4(3)]) > tolerance4)) &
        & call print_and_register( 'send strided real kind=4 to kind=4 self failed.')

      co_real_k8 = -1.0
      co_real_k8(::2)[2] = real_k8(1:3)
      print *, co_real_k8
      if (any(abs(co_real_k8 - [real_k8(1), REAL(-1.0, 8), real_k8(2), REAL(-1.0, 8), real_k8(3)]) > tolerance8)) &
        & call print_and_register( 'send strided real kind=8 to kind=8 self failed.')

      co_real_k8 = -2.0
      co_real_k8(2::2)[2] = real_k4(4:5)
      print *, co_real_k8(1:5), lbound(real_k8, 1)
      if (any(abs(co_real_k8 - [REAL(-2.0, 8), real_k8(4), REAL(-2.0, 8), real_k8(5), REAL(-2.0, 8)]) > tolerance4to8)) &
        & call print_and_register( 'send strided real kind=4 to kind=8 self failed.')

      co_real_k4 = -2.0
      co_real_k4(2::2)[2] = real_k8(0:1)
      print *, co_real_k4
      if (any(abs(co_real_k4 - [-2.0, real_k4(0), -2.0, real_k4(1), -2.0]) > tolerance4)) &
        & call print_and_register( 'send strided real kind=8 to kind=4 self failed.')

      co_real_k4(1:5) = real_k4(5:1:-1)
      co_real_k4(::2)[2] = co_real_k4(3:1:-1)
      print *, co_real_k4
      if (any(abs(co_real_k4 - [real_k4(3), real_k4(4), real_k4(4), real_k4(2), real_k4(5)]) > tolerance4)) &
        & call print_and_register( 'send strided with temp real kind=4 to kind=4 self failed.')

      co_real_k8(1:5) = real_k8(5:1:-1)
      co_real_k8(::2)[2] = co_real_k8(3:1:-1)
      print *, co_real_k8
      if (any(abs(co_real_k8 - [real_k8(3), real_k8(4), real_k8(4), real_k8(2), real_k8(5)]) > tolerance8)) &
        & call print_and_register( 'send strided with temp real kind=8 to kind=8 self failed.')
    end if

    ! Transfer to other image now.
    sync all
    co_int_k4 = -1
    co_int_k1 = INT(-1, 1)
    co_real_k8 = -1.0
    co_real_k4 = REAL(-1.0, 4)
    sync all
    if (me == 1) then
      co_int_k4(::2)[2] = [ 15, 13, 11]

      co_int_k1(::2)[2] = [INT(-15, 1), INT(-13, 1), INT(-11, 1)]

      co_real_k8(::2)[2] = [REAL(1.3, 8), REAL(1.5, 8), REAL(1.7, 8)]

      co_real_k4(::2)[2] = [REAL(1.3, 4), REAL(1.5, 4), REAL(1.7, 4)]
    end if

    sync all
    if (me == 2) then
      print *, co_int_k4
      if (any(co_int_k4 /= [15, -1, 13, -1, 11])) &
        call print_and_register( 'strided send int kind=4 to kind=4 to image 2 failed.')

      print *, co_int_k1
      if (any(co_int_k1 /= [INT(-15, 1), INT(-1, 1), INT(-13, 1), INT(-1, 1), INT(-11, 1)])) &
        & call print_and_register( 'strided send int kind=1 to kind=1 to image 2 failed.')

      print *, co_real_k8
      if (any(abs(co_real_k8 - [1.3, -1.0, 1.5, -1.0, 1.7]) > tolerance8)) &
        & call print_and_register( 'strided send real kind=8 to kind=8 to image 2 failed.')

      print *, co_real_k4
      if (any(abs(co_real_k4 - [REAL(1.3, 4), REAL(-1.0, 4), REAL(1.5, 4), REAL(-1.0, 4), REAL(1.7, 4)]) > tolerance4)) &
        & call print_and_register( 'strided send real kind=4 to kind=4 to image 2 failed.')
    end if

    ! now with strides and kind conversion
    sync all
    co_int_k4 = -1
    co_int_k1 = INT(-1, 1)
    co_real_k8 = -1.0
    co_real_k4 = REAL(-1.0, 4)
    sync all
    if (me == 1) then
      co_int_k4(::2)[2] = [INT(15, 1), INT(13, 1), INT(11, 1)]

      co_int_k1(::2)[2] = [-15, -13, -11]

      co_real_k8(::2)[2] = [REAL(1.3, 4), REAL(1.5, 4), REAL(1.7, 4)]

      co_real_k4(::2)[2] = [REAL(1.3, 8), REAL(1.5, 8), REAL(1.7, 8)]
    end if

    sync all
    if (me == 2) then
      print *, co_int_k4
      if (any(co_int_k4 /= [15, -1, 13, -1, 11])) &
        call print_and_register( 'strided send int kind=1 to kind=4 to image 2 failed.')

      print *, co_int_k1
      if (any(co_int_k1 /= [INT(-15, 1), INT(-1, 1), INT(-13, 1), INT(-1, 1), INT(-11, 1)])) &
        & call print_and_register( 'strided send int kind=4 to kind=1 to image 2 failed.')

      print *, co_real_k8
      if (any(abs(co_real_k8 - [1.3, -1.0, 1.5, -1.0, 1.7]) > tolerance8)) &
        & call print_and_register( 'strided send real kind=4 to kind=8 to image 2 failed.')

      print *, co_real_k4
      if (any(abs(co_real_k4 - [REAL(1.3, 4), REAL(-1.0, 4), REAL(1.5, 4), REAL(-1.0, 4), REAL(1.7, 4)]) > tolerance4)) &
        & call print_and_register( 'strided send real kind=8 to kind=4 to image 2 failed.')
    end if

    ! now with strides and type conversion
    sync all
    co_int_k4 = -1
    co_int_k1 = INT(-1, 1)
    co_real_k8 = -1.0
    co_real_k4 = REAL(-1.0, 4)
    sync all
    if (me == 1) then
      co_int_k4(::2)[2] = [15.0, 13.0, 11.0]

      co_int_k1(::2)[2] = [-15.0, -13.0, -11.0]

      co_real_k8(::2)[2] = [13, 15, 17]

      co_real_k4(::2)[2] = [23, 25, 27]
    end if

    sync all
    if (me == 2) then
      print *, co_int_k4
      if (any(co_int_k4 /= [15, -1, 13, -1, 11])) &
        call print_and_register( 'strided send real kind=4 to int kind=4 to image 2 failed.')

      print *, co_int_k1
      if (any(co_int_k1 /= [INT(-15, 1), INT(-1, 1), INT(-13, 1), INT(-1, 1), INT(-11, 1)])) &
        & call print_and_register( 'strided send int real kind=4 to int kind=1 to image 2 failed.')

      print *, co_real_k8
      if (any(abs(co_real_k8 - [13.0, -1.0, 15.0, -1.0, 17.0]) > tolerance8)) &
        & call print_and_register( 'strided send int kind=4 to real kind=8 to image 2 failed.')

      print *, co_real_k4
      if (any(abs(co_real_k4 - [REAL(23, 4), REAL(-1.0, 4), REAL(25, 4), REAL(-1.0, 4), REAL(27, 4)]) > tolerance4)) &
        & call print_and_register( 'strided send int kind=4 to real kind=4 to image 2 failed.')
    end if

    select case(me)
     case(1)
       if (error_printed) error stop
       sync images(2)
       print *, "Test passed."
     case(2)
       sync images(1)
       if (error_printed) error stop
    end select
    sync all
  end associate

contains

  subroutine print_and_register(error_message)
    use iso_fortran_env, only : error_unit
    character(len=*), intent(in) :: error_message
    write(error_unit,*) error_message
    error_printed=.true.
  end subroutine

end program send_convert_nums

! vim:ts=2:sts=2:sw=2:


text:
cc ------------ dpmjet3.4 - authors: S.Roesler, R.Engel, J.Ranft -------
cc -------- phojet1.12-40 - authors: S.Roesler, R.Engel, J.Ranft -------
cc                                                      - oct'13 -------
cc ----------- pythia-6.4 - authors: Torbjorn Sjostrand, Lund'10 -------
cc ---------------------------------------------------------------------
cc                                  converted for use with FLUKA -------
cc                                                      - oct'13 -------

C...PYSFDC
C...Calculates decays of sfermions.

      SUBROUTINE PYSFDC(KFIN,XLAM,IDLAM,IKNT)

C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)
      INTEGER PYCOMP
C...Parameter statement to help give large particle numbers.
      PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
     &KEXCIT=4000000,KDIMEN=5000000)
C...Commonblocks.
      include 'inc/pydat1'
      include 'inc/pydat2'
      include 'inc/pymssm'
      include 'inc/pyssmt'

C...Local variables.
      COMPLEX*16 ZMIXC(4,4),VMIXC(2,2),UMIXC(2,2)
      COMPLEX*16 CAL,CAR,CBL,CBR,CALP,CARP,CBLP,CBRP,CA,CB
      INTEGER KFIN,KCIN
      DOUBLE PRECISION XMI,XMJ,XMF,XMSF1,XMSF2,XMW,XMW2,XMZ,AXMJ
      DOUBLE PRECISION XMI2,XMI3,XMA2,XMB2,XMFP
      DOUBLE PRECISION PYLAMF,XL
      DOUBLE PRECISION TANW,XW,AEM,C1,AS
      DOUBLE PRECISION AL,AR,BL,BR
      DOUBLE PRECISION CH1,CH2,CH3,CH4
      DOUBLE PRECISION XMBOT,XMTOP
      DOUBLE PRECISION XLAM(0:400)
      INTEGER IDLAM(400,3)
      INTEGER LKNT,IX,ILR,IDU,J,I,IKNT,IFL,II
      DOUBLE PRECISION SR2
      DOUBLE PRECISION CBETA,SBETA
      DOUBLE PRECISION CW
      DOUBLE PRECISION BETA,ALFA,XMU,AT,AB,ATRIT,ATRIB,ATRIL
      DOUBLE PRECISION COSA,SINA,TANB
      DOUBLE PRECISION PYALEM,PI,PYALPS,EI
      DOUBLE PRECISION GHRR,GHLL,GHLR,XMB,BLR
      INTEGER IG,KF1,KF2
      INTEGER IGG(4),KFNCHI(4),KFCCHI(2)
      DATA IGG/23,25,35,36/
      DATA PI/3.141592653589793238462643383279D+00/
      DATA SR2/1.414213562373095048801688724210D+00/
      DATA KFNCHI/1000022,1000023,1000025,1000035/
      DATA KFCCHI/1000024,1000037/

C...COUNT THE NUMBER OF DECAY MODES
      LKNT=0

C...NO NU_R DECAYS
      IF(KFIN.EQ.KSUSY2+12.OR.KFIN.EQ.KSUSY2+14.OR.
     &KFIN.EQ.KSUSY2+16) RETURN

      XMW=PMAS(24,1)
      XMW2=XMW**2
      XMZ=PMAS(23,1)
      XW=PARU(102)
      TANW = SQRT(XW/(1D0-XW))
      CW=SQRT(1D0-XW)

      DO 110 I=1,4
        DO 100 J=1,4
          ZMIXC(J,I)=CMPLX(ZMIX(J,I),ZMIXI(J,I))
  100   CONTINUE
  110 CONTINUE
      DO 130 I=1,2
        DO 120 J=1,2
           VMIXC(J,I)=CMPLX(VMIX(J,I),VMIXI(J,I))
           UMIXC(J,I)=CMPLX(UMIX(J,I),UMIXI(J,I))
  120   CONTINUE
  130 CONTINUE

C...KCIN
      KCIN=PYCOMP(KFIN)
C...ILR is 1 for left and 2 for right.
      ILR=KFIN/KSUSY1
C...IFL is matching non-SUSY flavour.
      IFL=MOD(KFIN,KSUSY1)
C...IDU is weak isospin, 1 for down and 2 for up.
      IDU=2-MOD(IFL,2)

      XMI=PMAS(KCIN,1)
      XMI2=XMI**2
      AEM=PYALEM(XMI2)
      AS =PYALPS(XMI2)
      C1=AEM/XW
      XMI3=XMI**3
      EI=KCHG(IFL,1)/3D0

      XMBOT=PYMRUN(5,XMI2)
      XMTOP=PYMRUN(6,XMI2)

      TANB=RMSS(5)
      BETA=ATAN(TANB)
      ALFA=RMSS(18)
      CBETA=COS(BETA)
      SBETA=TANB*CBETA
      SINA=SIN(ALFA)
      COSA=COS(ALFA)
      XMU=-RMSS(4)
      ATRIT=RMSS(16)
      ATRIB=RMSS(15)
      ATRIL=RMSS(17)

C...2-BODY DECAYS OF SFERMION -> GRAVITINO + FERMION

      IF(IMSS(11).EQ.1) THEN
        XMP=RMSS(29)
        IDG=39+KSUSY1
        XMGR=PMAS(PYCOMP(IDG),1)
        XFAC=(XMI2/(XMP*XMGR))**2*XMI/48D0/PI
        IF(IFL.EQ.5) THEN
          XMF=XMBOT
        ELSEIF(IFL.EQ.6) THEN
          XMF=XMTOP
        ELSE
          XMF=PMAS(IFL,1)
        ENDIF
        IF(XMI.GT.XMGR+XMF) THEN
          LKNT=LKNT+1
          IDLAM(LKNT,1)=IDG
          IDLAM(LKNT,2)=IFL
          IDLAM(LKNT,3)=0
          XLAM(LKNT)=XFAC*(1D0-XMF**2/XMI2)**4
        ENDIF
      ENDIF

C...2-BODY DECAYS OF SFERMION -> FERMION + GAUGE/GAUGINO

C...CHARGED DECAYS:
      DO 140 IX=1,2
C...DI -> U CHI1-,CHI2-
        IF(IDU.EQ.1) THEN
          XMFP=PMAS(IFL+1,1)
          XMF =PMAS(IFL,1)
C...UI -> D CHI1+,CHI2+
        ELSE
          XMFP=PMAS(IFL-1,1)
          XMF =PMAS(IFL,1)
        ENDIF
        XMJ=SMW(IX)
        AXMJ=ABS(XMJ)
        IF(XMI.GE.AXMJ+XMFP) THEN
          XMA2=XMJ**2
          XMB2=XMFP**2
          IF(IDU.EQ.2) THEN
            IF(IFL.EQ.6) THEN
              XMFP=XMBOT
              XMF =XMTOP
            ELSEIF(IFL.LT.6) THEN
              XMF=0D0
              XMFP=0D0
            ENDIF
            CBL=VMIXC(IX,1)
            CAL=-XMFP*UMIXC(IX,2)/SR2/XMW/CBETA
            CBR=-XMF*VMIXC(IX,2)/SR2/XMW/SBETA
            CAR=0D0
          ELSE
            IF(IFL.EQ.5) THEN
              XMF =XMBOT
              XMFP=XMTOP
            ELSEIF(IFL.LT.5) THEN
              XMF=0D0
              XMFP=0D0
            ENDIF
            CBL=UMIXC(IX,1)
            CAL=-XMFP*VMIXC(IX,2)/SR2/XMW/SBETA
            CBR=-XMF*UMIXC(IX,2)/SR2/XMW/CBETA
            CAR=0D0
          ENDIF

          CALP=SFMIX(IFL,1)*CAL + SFMIX(IFL,2)*CAR
          CBLP=SFMIX(IFL,1)*CBL + SFMIX(IFL,2)*CBR
          CARP=SFMIX(IFL,4)*CAR + SFMIX(IFL,3)*CAL
          CBRP=SFMIX(IFL,4)*CBR + SFMIX(IFL,3)*CBL
          CAL=CALP
          CBL=CBLP
          CAR=CARP
          CBR=CBRP

C...F1 -> F` CHI
          IF(ILR.EQ.1) THEN
            CA=CAL
            CB=CBL
C...F2 -> F` CHI
          ELSE
            CA=CAR
            CB=CBR
          ENDIF
          LKNT=LKNT+1
          XL=PYLAMF(XMI2,XMA2,XMB2)
C...SPIN AVERAGE = 1/1 NOT 1/2....NO COLOR ENHANCEMENT
          XLAM(LKNT)=2D0*C1/8D0/XMI3*SQRT(XL)*((XMI2-XMB2-XMA2)*
     &    (ABS(CA)**2+ABS(CB)**2)-4D0*DBLE(CA*CONJG(CB))*XMJ*XMFP)
          IDLAM(LKNT,3)=0
          IF(IDU.EQ.1) THEN
            IDLAM(LKNT,1)=-KFCCHI(IX)
            IDLAM(LKNT,2)=IFL+1
          ELSE
            IDLAM(LKNT,1)=KFCCHI(IX)
            IDLAM(LKNT,2)=IFL-1
          ENDIF
        ENDIF
  140 CONTINUE

C...NEUTRAL DECAYS
      DO 150 IX=1,4
C...DI -> D CHI10
        XMF=PMAS(IFL,1)
        XMJ=SMZ(IX)
        AXMJ=ABS(XMJ)
        IF(XMI.GE.AXMJ+XMF) THEN
          XMA2=XMJ**2
          XMB2=XMF**2
          IF(IDU.EQ.1) THEN
            IF(IFL.EQ.5) THEN
              XMF=XMBOT
            ELSEIF(IFL.LT.5) THEN
              XMF=0D0
            ENDIF
            CBL=-ZMIXC(IX,2)+TANW*ZMIXC(IX,1)*(2D0*EI+1)
            CAL=XMF*ZMIXC(IX,3)/XMW/CBETA
            CAR=-2D0*EI*TANW*ZMIXC(IX,1)
            CBR=CAL
          ELSE
            IF(IFL.EQ.6) THEN
              XMF=XMTOP
            ELSEIF(IFL.LT.5) THEN
              XMF=0D0
            ENDIF
            CBL=ZMIXC(IX,2)+TANW*ZMIXC(IX,1)*(2D0*EI-1)
            CAL=XMF*ZMIXC(IX,4)/XMW/SBETA
            CAR=-2D0*EI*TANW*ZMIXC(IX,1)
            CBR=CAL
          ENDIF

          CALP=SFMIX(IFL,1)*CAL + SFMIX(IFL,2)*CAR
          CBLP=SFMIX(IFL,1)*CBL + SFMIX(IFL,2)*CBR
          CARP=SFMIX(IFL,4)*CAR + SFMIX(IFL,3)*CAL
          CBRP=SFMIX(IFL,4)*CBR + SFMIX(IFL,3)*CBL
          CAL=CALP
          CBL=CBLP
          CAR=CARP
          CBR=CBRP

C...F1 -> F CHI
          IF(ILR.EQ.1) THEN
            CA=CAL
            CB=CBL
C...F2 -> F CHI
          ELSE
            CA=CAR
            CB=CBR
          ENDIF
          LKNT=LKNT+1
          XL=PYLAMF(XMI2,XMA2,XMB2)
C...SPIN AVERAGE = 1/1 NOT 1/2....NO COLOR ENHANCEMENT
          XLAM(LKNT)=C1/8D0/XMI3*SQRT(XL)*((XMI2-XMB2-XMA2)*
     &    (ABS(CA)**2+ABS(CB)**2)-4D0*DBLE(CA*CONJG(CB))*XMJ*XMF)
          IDLAM(LKNT,1)=KFNCHI(IX)
          IDLAM(LKNT,2)=IFL
          IDLAM(LKNT,3)=0
        ENDIF
  150 CONTINUE

C...2-BODY DECAYS TO SM GAUGE AND HIGGS BOSONS
C...IG=23,25,35,36
      DO 160 II=1,4
        IG=IGG(II)
        IF(ILR.EQ.1) GOTO 160
        XMB=PMAS(IG,1)
        XMSF1=PMAS(PYCOMP(KFIN-KSUSY1),1)
        IF(XMI.LT.XMSF1+XMB) GOTO 160
        IF(IG.EQ.23) THEN
          BL=-SIGN(.5D0,EI)/CW+EI*XW/CW
          BR=EI*XW/CW
          BLR=0D0
        ELSEIF(IG.EQ.25) THEN
          IF(IFL.EQ.5) THEN
            XMF=XMBOT
          ELSEIF(IFL.EQ.6) THEN
            XMF=XMTOP
          ELSEIF(IFL.LT.5) THEN
            XMF=0D0
          ELSE
            XMF=PMAS(IFL,1)
          ENDIF
          IF(IDU.EQ.2) THEN
            GHLL=XMZ/CW*(0.5D0-EI*XW)*(-SIN(ALFA+BETA))+
     &      XMF**2/XMW*COSA/SBETA
            GHRR=XMZ/CW*(EI*XW)*(-SIN(ALFA+BETA))+
     &      XMF**2/XMW*COSA/SBETA
          ELSE
            GHLL=XMZ/CW*(0.5D0-EI*XW)*(-SIN(ALFA+BETA))+
     &      XMF**2/XMW*(-SINA)/CBETA
            GHRR=XMZ/CW*(EI*XW)*(-SIN(ALFA+BETA))+
     &      XMF**2/XMW*(-SINA)/CBETA
          ENDIF
          IF(IFL.EQ.5) THEN
            AT=ATRIB
          ELSEIF(IFL.EQ.6) THEN
            AT=ATRIT
          ELSEIF(IFL.EQ.15) THEN
            AT=ATRIL
          ELSE
            AT=0D0
          ENDIF
C.........need to complexify
          IF(IDU.EQ.2) THEN
            GHLR=XMF/2D0/XMW/SBETA*(-XMU*SINA+
     &      AT*COSA)
          ELSE
            GHLR=XMF/2D0/XMW/CBETA*(XMU*COSA-
     &      AT*SINA)
          ENDIF
          BL=GHLL
          BR=GHRR
          BLR=-GHLR
        ELSEIF(IG.EQ.35) THEN
          IF(IFL.EQ.5) THEN
            XMF=XMBOT
          ELSEIF(IFL.EQ.6) THEN
            XMF=XMTOP
          ELSEIF(IFL.LT.5) THEN
            XMF=0D0
          ELSE
            XMF=PMAS(IFL,1)
          ENDIF
          IF(IDU.EQ.2) THEN
            GHLL=XMZ/CW*(0.5D0-EI*XW)*COS(ALFA+BETA)+
     &      XMF**2/XMW*SINA/SBETA
            GHRR=XMZ/CW*(EI*XW)*COS(ALFA+BETA)+
     &      XMF**2/XMW*SINA/SBETA
          ELSE
            GHLL=XMZ/CW*(0.5D0-EI*XW)*COS(ALFA+BETA)+
     &      XMF**2/XMW*COSA/CBETA
            GHRR=XMZ/CW*(EI*XW)*COS(ALFA+BETA)+
     &      XMF**2/XMW*COSA/CBETA
          ENDIF
          IF(IFL.EQ.5) THEN
            AT=ATRIB
          ELSEIF(IFL.EQ.6) THEN
            AT=ATRIT
          ELSEIF(IFL.EQ.15) THEN
            AT=ATRIL
          ELSE
            AT=0D0
          ENDIF
C.........Need to complexify
          IF(IDU.EQ.2) THEN
            GHLR=XMF/2D0/XMW/SBETA*(XMU*COSA+
     &      AT*SINA)
          ELSE
            GHLR=XMF/2D0/XMW/CBETA*(XMU*SINA+
     &      AT*COSA)
          ENDIF
          BL=GHLL
          BR=GHRR
          BLR=GHLR
        ELSEIF(IG.EQ.36) THEN
          GHLL=0D0
          GHRR=0D0
          IF(IFL.EQ.5) THEN
            XMF=XMBOT
          ELSEIF(IFL.EQ.6) THEN
            XMF=XMTOP
          ELSEIF(IFL.LT.5) THEN
            XMF=0D0
          ELSE
            XMF=PMAS(IFL,1)
          ENDIF
          IF(IFL.EQ.5) THEN
            AT=ATRIB
          ELSEIF(IFL.EQ.6) THEN
            AT=ATRIT
          ELSEIF(IFL.EQ.15) THEN
            AT=ATRIL
          ELSE
            AT=0D0
          ENDIF
C.........Need to complexify
          IF(IDU.EQ.2) THEN
            GHLR=XMF/2D0/XMW*(-XMU+AT/TANB)
          ELSE
            GHLR=XMF/2D0/XMW/(-XMU+AT*TANB)
          ENDIF
          BL=GHLL
          BR=GHRR
          BLR=GHLR
        ENDIF
        AL=SFMIX(IFL,1)*SFMIX(IFL,3)*BL+
     &  SFMIX(IFL,2)*SFMIX(IFL,4)*BR+
     &  (SFMIX(IFL,1)*SFMIX(IFL,4)+SFMIX(IFL,3)*SFMIX(IFL,2))*BLR
        XL=PYLAMF(XMI2,XMSF1**2,XMB**2)
        LKNT=LKNT+1
        IF(IG.EQ.23) THEN
          XLAM(LKNT)=C1/4D0/XMI3*XL**1.5D0/XMB**2*AL**2
        ELSE
          XLAM(LKNT)=C1/4D0/XMI3*SQRT(XL)*AL**2
        ENDIF
        IDLAM(LKNT,3)=0
        IDLAM(LKNT,1)=KFIN-KSUSY1
        IDLAM(LKNT,2)=IG
  160 CONTINUE

C...SF -> SF' + W
      XMB=PMAS(24,1)
      IF(MOD(IFL,2).EQ.0) THEN
        KF1=KSUSY1+IFL-1
      ELSE
        KF1=KSUSY1+IFL+1
      ENDIF
      KF2=KF1+KSUSY1
      XMSF1=PMAS(PYCOMP(KF1),1)
      XMSF2=PMAS(PYCOMP(KF2),1)
      IF(XMI.GT.XMB+XMSF1) THEN
        IF(MOD(IFL,2).EQ.0) THEN
          IF(ILR.EQ.1) THEN
            AL=1D0/SR2*SFMIX(IFL,1)*SFMIX(IFL-1,1)
          ELSE
            AL=1D0/SR2*SFMIX(IFL,3)*SFMIX(IFL-1,1)
          ENDIF
        ELSE
          IF(ILR.EQ.1) THEN
            AL=1D0/SR2*SFMIX(IFL,1)*SFMIX(IFL+1,1)
          ELSE
            AL=1D0/SR2*SFMIX(IFL,3)*SFMIX(IFL+1,1)
          ENDIF
        ENDIF
        XL=PYLAMF(XMI2,XMSF1**2,XMB**2)
        LKNT=LKNT+1
        XLAM(LKNT)=C1/4D0/XMI3*XL**1.5D0/XMB**2*AL**2
        IDLAM(LKNT,3)=0
        IDLAM(LKNT,1)=KF1
        IDLAM(LKNT,2)=SIGN(24,KCHG(IFL,1))
      ENDIF
      IF(XMI.GT.XMB+XMSF2) THEN
        IF(MOD(IFL,2).EQ.0) THEN
          IF(ILR.EQ.1) THEN
            AL=1D0/SR2*SFMIX(IFL,1)*SFMIX(IFL-1,3)
          ELSE
            AL=1D0/SR2*SFMIX(IFL,3)*SFMIX(IFL-1,3)
          ENDIF
        ELSE
          IF(ILR.EQ.1) THEN
            AL=1D0/SR2*SFMIX(IFL,1)*SFMIX(IFL+1,3)
          ELSE
            AL=1D0/SR2*SFMIX(IFL,3)*SFMIX(IFL+1,3)
          ENDIF
        ENDIF
        XL=PYLAMF(XMI2,XMSF2**2,XMB**2)
        LKNT=LKNT+1
        XLAM(LKNT)=C1/4D0/XMI3*XL**1.5D0/XMB**2*AL**2
        IDLAM(LKNT,3)=0
        IDLAM(LKNT,1)=KF2
        IDLAM(LKNT,2)=SIGN(24,KCHG(IFL,1))
      ENDIF

C...SF -> SF' + HC
      XMB=PMAS(37,1)
      IF(MOD(IFL,2).EQ.0) THEN
        KF1=KSUSY1+IFL-1
      ELSE
        KF1=KSUSY1+IFL+1
      ENDIF
      KF2=KF1+KSUSY1
      XMSF1=PMAS(PYCOMP(KF1),1)
      XMSF2=PMAS(PYCOMP(KF2),1)
      IF(XMI.GT.XMB+XMSF1) THEN
        XMF=0D0
        XMFP=0D0
        AT=0D0
        AB=0D0
        IF(MOD(IFL,2).EQ.0) THEN
C...T1-> B1 HC
          IF(ILR.EQ.1) THEN
            CH1=-SFMIX(IFL,1)*SFMIX(IFL-1,1)
            CH2= SFMIX(IFL,2)*SFMIX(IFL-1,2)
            CH3=-SFMIX(IFL,1)*SFMIX(IFL-1,2)
            CH4=-SFMIX(IFL,2)*SFMIX(IFL-1,1)
C...T2-> B1 HC
          ELSE
            CH1= SFMIX(IFL,3)*SFMIX(IFL-1,1)
            CH2=-SFMIX(IFL,4)*SFMIX(IFL-1,2)
            CH3= SFMIX(IFL,3)*SFMIX(IFL-1,2)
            CH4= SFMIX(IFL,4)*SFMIX(IFL-1,1)
          ENDIF
          IF(IFL.EQ.6) THEN
            XMF=XMTOP
            XMFP=XMBOT
            AT=ATRIT
            AB=ATRIB
          ENDIF
        ELSE
C...B1 -> T1 HC
          IF(ILR.EQ.1) THEN
            CH1=-SFMIX(IFL+1,1)*SFMIX(IFL,1)
            CH2= SFMIX(IFL+1,2)*SFMIX(IFL,2)
            CH3=-SFMIX(IFL+1,1)*SFMIX(IFL,2)
            CH4=-SFMIX(IFL+1,2)*SFMIX(IFL,1)
C...B2-> T1 HC
          ELSE
            CH1= SFMIX(IFL,3)*SFMIX(IFL+1,1)
            CH2=-SFMIX(IFL,4)*SFMIX(IFL+1,2)
            CH3= SFMIX(IFL,4)*SFMIX(IFL+1,1)
            CH4= SFMIX(IFL,3)*SFMIX(IFL+1,2)
          ENDIF
          IF(IFL.EQ.5) THEN
            XMF=XMTOP
            XMFP=XMBOT
            AT=ATRIT
            AB=ATRIB
          ENDIF
        ENDIF
        XL=PYLAMF(XMI2,XMSF1**2,XMB**2)
        LKNT=LKNT+1
C.......Need to complexify
        AL=CH1*(XMW2*2D0*CBETA*SBETA-XMFP**2*TANB-XMF**2/TANB)+
     &  CH2*2D0*XMF*XMFP/(2D0*CBETA*SBETA)+
     &  CH3*XMFP*(-XMU+AB*TANB)+CH4*XMF*(-XMU+AT/TANB)
        XLAM(LKNT)=C1/8D0/XMI3*SQRT(XL)/XMW2*AL**2
        IDLAM(LKNT,3)=0
        IDLAM(LKNT,1)=KF1
        IDLAM(LKNT,2)=SIGN(37,KCHG(IFL,1))
      ENDIF
      IF(XMI.GT.XMB+XMSF2) THEN
        XMF=0D0
        XMFP=0D0
        AT=0D0
        AB=0D0
        IF(MOD(IFL,2).EQ.0) THEN
C...T1-> B2 HC
          IF(ILR.EQ.1) THEN
            CH1= SFMIX(IFL-1,3)*SFMIX(IFL,1)
            CH2=-SFMIX(IFL-1,4)*SFMIX(IFL,2)
            CH3= SFMIX(IFL-1,4)*SFMIX(IFL,1)
            CH4= SFMIX(IFL-1,3)*SFMIX(IFL,2)
C...T2-> B2 HC
          ELSE
            CH1= -SFMIX(IFL,3)*SFMIX(IFL-1,3)
            CH2= SFMIX(IFL,4)*SFMIX(IFL-1,4)
            CH3= -SFMIX(IFL,3)*SFMIX(IFL-1,4)
            CH4= -SFMIX(IFL,4)*SFMIX(IFL-1,3)
          ENDIF
          IF(IFL.EQ.6) THEN
            XMF=XMTOP
            XMFP=XMBOT
            AT=ATRIT
            AB=ATRIB
          ENDIF
        ELSE
C...B1 -> T2 HC
          IF(ILR.EQ.1) THEN
            CH1= SFMIX(IFL+1,3)*SFMIX(IFL,1)
            CH2=-SFMIX(IFL+1,4)*SFMIX(IFL,2)
            CH3= SFMIX(IFL+1,3)*SFMIX(IFL,2)
            CH4= SFMIX(IFL+1,4)*SFMIX(IFL,1)
C...B2-> T2 HC
          ELSE
            CH1= -SFMIX(IFL+1,3)*SFMIX(IFL,3)
            CH2= SFMIX(IFL+1,4)*SFMIX(IFL,4)
            CH3= -SFMIX(IFL+1,3)*SFMIX(IFL,4)
            CH4= -SFMIX(IFL+1,4)*SFMIX(IFL,3)
          ENDIF
          IF(IFL.EQ.5) THEN
            XMF=XMTOP
            XMFP=XMBOT
            AT=ATRIT
            AB=ATRIB
          ENDIF
        ENDIF
        XL=PYLAMF(XMI2,XMSF1**2,XMB**2)
        LKNT=LKNT+1
C.......Need to complexify
        AL=CH1*(XMW2*2D0*CBETA*SBETA-XMFP**2*TANB-XMF**2/TANB)+
     &  CH2*2D0*XMF*XMFP/(2D0*CBETA*SBETA)+
     &  CH3*XMFP*(-XMU+AB*TANB)+CH4*XMF*(-XMU+AT/TANB)
        XLAM(LKNT)=C1/8D0/XMI3*SQRT(XL)/XMW2*AL**2
        IDLAM(LKNT,3)=0
        IDLAM(LKNT,1)=KF2
        IDLAM(LKNT,2)=SIGN(37,KCHG(IFL,1))
      ENDIF

C...2-BODY DECAYS OF SQUARK -> QUARK GLUINO

      IF(IFL.LE.6) THEN
        XMFP=0D0
        XMF=0D0
        IF(IFL.EQ.6) XMF=PMAS(6,1)
        IF(IFL.EQ.5) XMF=PMAS(5,1)
        XMJ=PMAS(PYCOMP(KSUSY1+21),1)
        AXMJ=ABS(XMJ)
        IF(XMI.GE.AXMJ+XMF) THEN
          AL=-SFMIX(IFL,3)
          BL=SFMIX(IFL,1)
          AR=-SFMIX(IFL,4)
          BR=SFMIX(IFL,2)
C...F1 -> F CHI
          IF(ILR.EQ.1) THEN
            XCA=AL
            XCB=BL
C...F2 -> F CHI
          ELSE
            XCA=AR
            XCB=BR
          ENDIF
          LKNT=LKNT+1
          XMA2=XMJ**2
          XMB2=XMF**2
          XL=PYLAMF(XMI2,XMA2,XMB2)
          XLAM(LKNT)=4D0/3D0*AS/2D0/XMI3*SQRT(XL)*((XMI2-XMB2-XMA2)*
     &    (XCA**2+XCB**2)+4D0*XCA*XCB*XMJ*XMF)
          IDLAM(LKNT,1)=KSUSY1+21
          IDLAM(LKNT,2)=IFL
          IDLAM(LKNT,3)=0
        ENDIF
      ENDIF

C...IF NOTHING ELSE FOR T1, THEN T1* -> C+CHI0
      IF(KFIN.EQ.KSUSY1+6.AND.PMAS(KCIN,1).GT.
     &PMAS(PYCOMP(KSUSY1+22),1)+PMAS(4,1)) THEN
C...THIS IS A BACK-OF-THE-ENVELOPE ESTIMATE
C...M = 1/(16PI**2)G**3 = G*2/(4PI) G/(4PI) = C1 * G/(4PI)
C...M*M = C1**2 * G**2/(16PI**2)
C...G = 1/(8PI)P/MI**2 * M*M = C1**3/(32PI**2)*LAM/(2*MI**3)
        LKNT=LKNT+1
        XL=PYLAMF(XMI2,0D0,PMAS(PYCOMP(KSUSY1+22),1)**2)
        XLAM(LKNT)=C1**3/64D0/PI**2/XMI3*SQRT(XL)
        IF(XLAM(LKNT).EQ.0) XLAM(LKNT)=1D-3
        IDLAM(LKNT,1)=KSUSY1+22
        IDLAM(LKNT,2)=4
        IDLAM(LKNT,3)=0
      ENDIF

C...R-violating sfermion decays (SKANDS).
      CALL PYRVSF(KFIN,XLAM,IDLAM,LKNT)

      IKNT=LKNT
      XLAM(0)=0D0
      DO 170 I=1,IKNT
        IF(XLAM(I).LT.0D0) XLAM(I)=0D0
        XLAM(0)=XLAM(0)+XLAM(I)
  170 CONTINUE
      IF(XLAM(0).EQ.0D0) XLAM(0)=1D-3

      RETURN
      END

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
          program demo_close
          implicit none
          character(len=256) :: message
          integer            :: ios
             open (10, file='employee.names', action='read', iostat=ios,iomsg=message)
             if (ios < 0) then
                ! perform end-of-file processing on the file connected to unit 10.

                close (10, status='keep',iostat=ios,iomsg=message)
                if(ios.ne.0)then
                   write(*,'(*(a))')'*demo_close* close error: ',trim(message)
                   stop 1
                endif
             elseif (ios > 0) then
                ! perform error processing on open
                write(*,'(*(a))')'*demo_close* open error: ',trim(message)
                stop 2
             endif
          end program demo_close

max_stars_repo_name: urbanjost/fortran-intrinsic-manpages

max_issues_repo_name: urbanjost/fortran-intrinsic-manpages

max_forks_repo_name: urbanjost/fortran-intrinsic-manpages