text:
	SUBROUTINE readinputfile(Nx,Ny,Nz,Nt,plotgap,Lx,Ly,Lz, &
								Es,DT,starttime,myid,ierr)
	!--------------------------------------------------------------------------------
	!
	!
	!  PURPOSE
	!
	!	Read inputfile intialize parameters, which are stocked in the Input File
	!
	! .. INPUT ..
	!	Nx				= number of modes in the x direction
	!	Ny				= number of modes in the y direction
	!	Nz				= number of modes in the z direction
	!	Nt				= the number of timesteps
	!	plotgap			= the number of timesteps to take before plotting
	!	myid			= number of MPI process
	!	ierr			= MPI error output variable
	!	Lx				= size of the periodic domain of computation in x direction
	!	Ly				= size of the periodic domain of computation in y direction
	!	Lz				= size of the periodic domain of computation in z direction
	!	DT				= the time step
	!	starttime		= initial time of computation
	!	InputFileName	= name of the Input File
	! REFERENCES
	!
	! ACCURACY
	!
	! ERROR INDICATORS AND WARNINGS
	!
	! FURTHER COMMENTS
	!---------------------------------------------------------------------------------
	! EXTERNAL ROUTINES REQUIRED
	IMPLICIT NONE
	INCLUDE 'mpif.h'
	! .. Scalar Arguments ..
	INTEGER(KIND=4), INTENT(IN)		::  myid
	INTEGER(KIND=4), INTENT(OUT)	::  Nx,Ny,Nz,Nt
	INTEGER(KIND=4), INTENT(OUT)	::  plotgap, ierr
	REAL(KIND=8), INTENT(OUT)		::  Lx, Ly, Lz, DT, starttime, Es
	! .. Local scalars ..
	INTEGER(KIND=4)					::  stat
	! .. Local Arrays ..
	CHARACTER*40					::  InputFileName
	INTEGER(KIND=4), DIMENSION(1:5)	::  intcomm
	REAL(KIND=8), DIMENSION(1:6)	::  dpcomm

	IF(myid.eq.0) THEN
		CALL GET_ENVIRONMENT_VARIABLE(NAME="inputfile",VALUE=InputFileName, STATUS=stat)
		IF(stat.NE.0) THEN
			PRINT*,"Set environment variable inputfile to the name of the"
			PRINT*,"file where the simulation parameters are set"
			STOP
		END IF
		OPEN(unit=11,FILE=trim(InputFileName),status="OLD")
		REWIND(11)
		READ(11,*) intcomm(1), intcomm(2), intcomm(3), intcomm(4), intcomm(5), &
			 	dpcomm(1), dpcomm(2), dpcomm(3), dpcomm(4), dpcomm(5), dpcomm(6)
		CLOSE(11)
		PRINT *,"NX ",intcomm(1)
		PRINT *,"NY ",intcomm(2)
		PRINT *,"NZ ",intcomm(3)
		PRINT *,"NT ",intcomm(4)
		PRINT *,"plotgap ",intcomm(5)
		PRINT *,"Lx ",dpcomm(1)
		PRINT *,"Ly ",dpcomm(2)
		PRINT *,"Lz ",dpcomm(3)
		PRINT *,"Es ",dpcomm(4)
		PRINT *,"Dt ",dpcomm(5)
		PRINT *,"strart time ",dpcomm(6)
		PRINT *,"Read inputfile"
	END IF
	CALL MPI_BCAST(dpcomm,6,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
	CALL MPI_BCAST(intcomm,5,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)

	Nx=intcomm(1)
	Ny=intcomm(2)
	Nz=intcomm(3)
	Nt=intcomm(4)
	plotgap=intcomm(5)
	Lx=dpcomm(1)
	Ly=dpcomm(2)
	Lz=dpcomm(3)
	Es=dpcomm(4)
	DT=dpcomm(5)
	starttime=dpcomm(6)

	END SUBROUTINE readinputfile

max_stars_repo_name: bcloutier/PSNM

max_issues_repo_name: bcloutier/PSNM

max_forks_repo_name: bcloutier/PSNM


        !COMPILER-GENERATED INTERFACE MODULE: Thu Oct 25 10:38:56 2018
        MODULE DECSF__genmod
          INTERFACE
            SUBROUTINE DECSF(X,Y,ZZ,FI0,TET0,FI,TET,H)
              REAL(KIND=4) :: X
              REAL(KIND=4) :: Y
              REAL(KIND=4) :: ZZ
              REAL(KIND=4) :: FI0
              REAL(KIND=4) :: TET0
              REAL(KIND=4) :: FI
              REAL(KIND=4) :: TET
              REAL(KIND=4) :: H
            END SUBROUTINE DECSF
          END INTERFACE
        END MODULE DECSF__genmod


text:
MODULE SOI_delta_scaling

	USE Type_Kinds, only : fpk => fp_kind

	IMPLICIT NONE

	REAL(fpk), PRIVATE, PARAMETER	:: ONE = 1.0_fpk
	REAL(fpk), PRIVATE, PARAMETER	:: TWO = 2.0_fpk
	REAL(fpk), PRIVATE, PARAMETER	:: ZERO = 0.0_fpk

PUBLIC  :: deltascaling, deltascaling_TL, deltascaling_AD

CONTAINS

	SUBROUTINE deltascaling( opd,  &  ! input
                              ssa,  &  ! input
                              asy,  &  ! input
                                f,  &  ! input
                             opd2,  &  ! output
                             ssa2,  &  ! output
                             asy2)     ! output

		REAL(fpk), intent(in), dimension(:)		::	opd, ssa, asy, f
		REAL(fpk), intent(out), dimension(:)	::	opd2, ssa2, asy2

		opd2 = opd * (ONE - ssa * f)
		ssa2 = (ONE-f)*ssa/(ONE - f*ssa)
		asy2 = (asy-f)/(ONE-f)

	END SUBROUTINE deltascaling

	SUBROUTINE deltascaling_TL( opd,   &  ! input
                                ssa,    &  ! input
                                asy,    &  ! input
                                  f,    &  ! input
						       opd_TL,  &  ! input
	                           ssa_TL,  &  ! input
                               asy_TL,  &  ! input
                                 f_TL,  &  ! input
                               opd2_TL, &  ! output
                               ssa2_TL, &  ! output
                               asy2_TL  )  ! output

		REAL(fpk), intent(in), dimension(:)		::	opd, ssa, asy, f
		REAL(fpk), intent(in), dimension(:)		::	opd_TL, ssa_TL, asy_TL, f_TL
		REAL(fpk), intent(out), dimension(:)	::	opd2_TL, ssa2_TL, asy2_TL

		opd2_TL = (ONE-ssa*f)*opd_TL - opd*(ssa*f_TL + f*ssa_TL)
		ssa2_TL = (ssa*(ssa-ONE)*f_TL + (ONE-f)*ssa_TL)/(ONE-ssa*f)**2
		asy2_TL = asy_TL / (ONE-f) - (ONE-asy)*f_TL/(ONE-f)**2

	END SUBROUTINE deltascaling_TL

	SUBROUTINE deltascaling_AD( opd,   &  ! input
                                ssa,    &  ! input
                                asy,    &  ! input
                                  f,    &  ! input
						       opd_AD,  &  ! inout
	                           ssa_AD,  &  ! inout
                               asy_AD,  &  ! inout
                                 f_AD,  &  ! inout
                               opd2_AD, &  ! inout
                               ssa2_AD, &  ! inout
                               asy2_AD  )  ! inout

		REAL(fpk), intent(in), dimension(:)     ::	opd, ssa, asy, f
		REAL(fpk), intent(inout), dimension(:)  ::	opd_AD, ssa_AD, asy_AD, f_AD
		REAL(fpk), intent(inout), dimension(:)  ::	opd2_AD, ssa2_AD, asy2_AD
		REAL(fpk), dimension(size(opd))         ::	denom1, denom2

		! useful constants
		denom1 = (ONE-ssa*f)
		denom1 = ONE/(denom1*denom1)
		denom2 = ONE/(ONE-f)

		! adjoint statements
		opd_AD = opd_AD + (ONE-ssa*f)*opd2_AD
		ssa_AD = ssa_AD + (ONE-f)*denom1 * ssa2_AD - f*opd*opd2_AD
		asy_AD = asy_AD + asy2_AD*denom2
		f_AD = f_AD - (ONE-asy)*denom2*denom2*asy2_AD + ssa*(ssa-ONE)*denom1*ssa2_AD - ssa*opd*opd2_AD
		opd2_AD = ZERO
		ssa2_AD = ZERO
		asy2_AD = ZERO

	END SUBROUTINE deltascaling_AD

END MODULE SOI_delta_scaling

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
      function diskdf5intez(vt,r,z)
      parameter (pi=3.1415926535)

      real psir0,psirz

      psir0=pot(r,0.0)
      if (z.eq.0.) then
         psirz=psir0
      else
         psirz=pot(r,z)
         endif
      ep=0.5*(vt*vt)-psir0
      am=r*vt
      ez=-psirz+psir0

      diskdf5intez=diskdf3intez(ep,am,ez)
      if(diskdf5intez.lt.0.) diskdf5intez = 0.

      return
      end

ext: f

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


!> @brief   Functions and types for string manipulation
module fplus_strings
    use, intrinsic :: iso_c_binding
    use fplus_object
    use fplus_hashcode
    implicit none

    !> @brief   A simple container for character strings
    type, extends(object) :: string
        !> @brief   the actual string content
        character (len=1000), public :: chars
    contains
        !> @brief   copy a c-string to this string object
        procedure, public :: set_cstr => string_set_cstr
        !> @brief   Returns a string representation of this object
        procedure, public :: to_string => string_to_string
        !> @brief   Calculate the hash code for this object
        procedure, public :: hashcode => string_hashcode
    end type

contains

    !> @public
    !> @brief       Calculate the hash code for this object
    function string_hashcode(this) result (res)
        class (string) :: this
        integer (kind=8) :: res

        ! local variables
        character (len=:), allocatable :: str
        integer (kind=1), dimension(:), allocatable :: chars
        integer :: l

        str = "String" // trim(this%chars)

        l = len(str)
        allocate(chars(l))
        chars = transfer(str, chars)
        call C_sdbm(chars, l, res)
    end function

    !> @public
    !> @brief       Returns a string representation of this object
    function string_to_string(this) result (res)
        class(string) :: this
        character (len=:), allocatable :: res
        res = trim(this%chars)
    end function

    !> @brief   determine the number of digits in a number, a negative number has one digit more
    function ndigits_of_integer(number) result(ndigits)
        integer :: ndigits
        class(*) :: number
        select type (number)
            type is (integer (kind=4))
                if (number > 0) then
                    ndigits = floor(log10(real(number)))+1
                else if (number == 0) then
                    ndigits = 1
                else if (number < 0) then
                    ndigits = floor(log10(real(abs(number))))+2
                end if
            type is (integer (kind=8))
                if (number > 0) then
                    ndigits = floor(log10(real(number)))+1
                else if (number == 0) then
                    ndigits = 1
                else if (number < 0) then
                    ndigits = floor(log10(real(abs(number))))+2
                end if
            class default
                ndigits = 0
        end select
    end function

    !> @brief   Convert all number typs to strings
    function number_to_string(number, format) result(res)
        character (len=:), allocatable :: res
        character (len=*), intent(in), optional :: format
        class(*) :: number

        ! local variables
        character (len=100) :: temp
        integer :: ndigits

        ! create the number string depending on the format
        select type (number)
            type is (real (kind=4))
                if (present(format)) then
                    write(temp, format) number
                else
                    write(temp, *) number
                end if
            type is (real (kind=8))
                if (present(format)) then
                    write(temp, format) number
                else
                    write(temp, *) number
                end if
            type is (integer (kind=4))
                if (present(format)) then
                    write(temp, format) number
                else
                    write(temp, *) number
                end if
            type is (integer (kind=8))
                if (present(format)) then
                    write(temp, format) number
                else
                    write(temp, *) number
                end if
            class default
                res = "unknown type in number_to_string"
        end select
        res = trim(adjustl(temp))
    end function

    !> @brief   Returns a string representation of intrinsic types and types that extend
    !>          the abstract type object with the to_string method
    function type_to_string(val) result(res)
        character (len=:), allocatable :: res
        class(*) :: val

        select type (val)
            type is (real (kind=4))
                res = number_to_string(val)
            type is (real (kind=8))
                res = number_to_string(val)
            type is (integer (kind=4))
                res = number_to_string(val)
            type is (integer (kind=8))
                res = number_to_string(val)
            type is (character (len=*))
                res = val
            type is (logical)
                if (val .eqv. .true.) then
                    res = "T"
                else
                    res = "F"
                end if
            class is (object)
                res = val%to_string()
            class default
                res = "unknown type in type_to_string"
        end select
    end function

    !> @brief   copy a c-string to this string object
    subroutine string_set_cstr(this, cstr)
        class(string) :: this
        character (kind=C_char) :: cstr(*)

        ! local variables
        integer :: i
        logical :: in_cstr

        ! loop over the length of this string
        in_cstr = .true.
        do i = 1, len(this%chars)
            if (cstr(i) == C_null_char) then
                in_cstr = .false.
            end if
            if (in_cstr) then
                this%chars(i:i) = cstr(i)
            else
                this%chars(i:i) = ' '
            end if
        end do
    end subroutine

end module fplus_strings


function gaussian_kernel(a, b) result(kernel)

    double precision, dimension(:), intent(in) :: a
    double precision, dimension(:), intent(in) :: b

    double precision :: kernel
    integer :: i
    double precision :: temp
    double precision :: sigma

    kernel = 0.0d0
    sigma = 724.0d0

    do i = 1, size(a)
        temp = a(i) - b(i)
        kernel = kernel + temp * temp
    enddo

    kernel = exp(-0.5d0 * sqrt(kernel) / (sigma*sigma))

end function gaussian_kernel

function laplace_kernel(a, b) result(kernel)

    double precision, dimension(:), intent(in) :: a
    double precision, dimension(:), intent(in) :: b

    double precision :: kernel

    kernel = sum(abs(a(:) - b(:)))
    kernel = exp(-0.00025118864315095795d0 * kernel)

end function laplace_kernel

subroutine kgaussian_kernel(a, na, b, nb, k, sigma)

    double precision, dimension(:,:), intent(in) :: a
    double precision, dimension(:,:), intent(in) :: b

    integer, intent(in) :: na, nb

    double precision, dimension(:,:), intent(inout) :: k
    double precision, intent(in) :: sigma

    double precision, allocatable, dimension(:) :: temp

    double precision :: inv_sigma
    integer :: l

    inv_sigma = -0.5d0 / (sigma*sigma)

    allocate(temp(size(a, dim=1)))

!$OMP PARALLEL DO PRIVATE(temp)
    do i = 1, nb
        do j = 1, na
            temp(:) = a(:,j) - b(:,i)
            K(j,i) = exp(inv_sigma * sqrt(sum(temp*temp)))
        enddo
    enddo
!$OMP END PARALLEL DO

    deallocate(temp)
end subroutine kgaussian_kernel

subroutine klaplace_kernel(a, na, b, nb, k, sigma)

    double precision, dimension(:,:), intent(in) :: a
    double precision, dimension(:,:), intent(in) :: b

    integer, intent(in) :: na, nb

    double precision, dimension(:,:), intent(inout) :: k
    double precision, intent(in) :: sigma

    double precision :: inv_sigma

    inv_sigma = -1.0d0 / sigma

!$OMP PARALLEL DO
    do i = 1, nb
        do j = 1, na
        k(j,i) = exp(inv_sigma * sum(abs(a(:,j) - b(:,i))))
        enddo
    enddo
!$OMP END PARALLEL DO

end subroutine klaplace_kernel


text:
      DOUBLE PRECISION FUNCTION DGENCHI(DF)
C**********************************************************************
C
C       DOUBLE PRECISION FUNCTION DGENCHI( DF )
C                Generate random value of CHIsquare variable
C
C
C                              Function
C
C
C     Generates random deviate from the distribution of a chisquare
C     with DF degrees of freedom random variable.
C
C
C                              Arguments
C
C
C     DF --> Degrees of freedom of the chisquare
C            (Must be positive)
C                    DOUBLE PRECISION DF
C
C
C                              Method
C
C
C     Uses relation between chisquare and gamma.
C
C**********************************************************************
C     .. Scalar Arguments ..
      DOUBLE PRECISION DF
C     ..
C     .. External Functions ..
      DOUBLE PRECISION DGENGAM
      EXTERNAL DGENGAM
C     ..
C     .. Executable Statements ..
      IF (.NOT. (DF.LE.0.0D0)) GO TO 10
      WRITE (*,FMT=*) 'DF <= 0 in DGENCHI - ABORT'
      WRITE (*,FMT=*) 'Value of DF: ',DF
      STOP 'DF <= 0 in DGENCHI - ABORT'

   10 DGENCHI = 2.0D0*DGENGAM(1.0D0,DF/2.0D0)
      RETURN

      END
c -------------------------------------------------
      DOUBLE PRECISION FUNCTION DGENGAM(A,R)
C**********************************************************************
C
C        DOUBLE PRECISION FUNCTION DGENGAM( A, R )
C           GENerates random deviates from GAMma distribution
C
C
C                              Function
C
C
C     Generates random deviates from the gamma distribution whose
C     density is
C          (A**R)/Gamma(R) * X**(R-1) * Exp(-A*X)
C
C
C                              Arguments
C
C
C     A --> Location parameter of Gamma distribution
C                         DOUBLE PRECISION A
C
C     R --> Shape parameter of Gamma distribution
C                         DOUBLE PRECISION R
C
C
C                              Method
C
C
C     Renames DSGAMMA from TOMS as slightly modified by BWB to use
C     DRANFD instead of SUNIF.
C
C     For details see:
C               (Case R >= 1.0)
C               Ahrens, J.H. and Dieter, U.
C               Generating Gamma Variates by a
C               Modified Rejection Technique.
C               Comm. ACM, 25,1 (Jan. 1982), 47 - 54.
C     Algorithm GD
C
C               (Case 0.0 <= R <= 1.0)
C               Ahrens, J.H. and Dieter, U.
C               Computer Methods for Sampling from Gamma,
C               Beta, Poisson and Binomial Distributions.
C               Computing, 12 (1974), 223-246/
C     Adapted algorithm GS.
C
C**********************************************************************
C     .. Scalar Arguments ..
      DOUBLE PRECISION A,R
C     ..
C     .. External Functions ..
      DOUBLE PRECISION DSGAMMA
      EXTERNAL DSGAMMA
C     ..
C     .. Executable Statements ..
      DGENGAM = DSGAMMA(R)
      DGENGAM = DGENGAM/A
      RETURN

      END
c -------------------------------------------------
      DOUBLE PRECISION FUNCTION DGENNOR(AV,SD)
C**********************************************************************
C
C        DOUBLE PRECISION FUNCTION DGENNOR( AV, SD )
C
C         GENerate random deviate from a NORmal distribution
C
C
C                              Function
C
C
C     Generates a single random deviate from a normal distribution
C     with mean, AV, and standard deviation, SD.
C
C
C                              Arguments
C
C
C     AV --> Mean of the normal distribution.
C                         DOUBLE PRECISION AV
C
C     SD --> Standard deviation of the normal distribution.
C                         DOUBLE PRECISION SD
C
C     DGENNOR <-- Generated normal deviate.
C                         DOUBLE PRECISION DGENNOR
C
C
C                              Method
C
C
C     Renames DSNORM from TOMS as slightly modified by BWB to use
C     DRANFD instead of SUNIF.
C
C     For details see:
C               Ahrens, J.H. and Dieter, U.
C               Extensions of Forsythe's Method for Random
C               Sampling from the Normal Distribution.
C               Math. Comput., 27,124 (Oct. 1973), 927 - 937.
C
C
C**********************************************************************
C     .. Scalar Arguments ..
      DOUBLE PRECISION AV,SD
C     ..
C     .. External Functions ..
      DOUBLE PRECISION DSNORM
      EXTERNAL DSNORM
C     ..
C     .. Executable Statements ..
      DGENNOR = SD*DSNORM() + AV
      RETURN

      END
c -------------------------------------------------
      DOUBLE PRECISION FUNCTION DGENUNF(LOW,HIGH)
C**********************************************************************
C
C         DOUBLE PRECISION FUNCTION DGENUNF( LOW, HIGH )
C
C               GeNerate Uniform Real between LOW and HIGH
C
C
C                              Function
C
C
C     Generates a real uniformly distributed between LOW and HIGH.
C
C
C                              Arguments
C
C
C     LOW --> Low bound (exclusive) on real value to be generated
C                    DOUBLE PRECISION LOW
C
C     HIGH --> High bound (exclusive) on real value to be generated
C                    DOUBLE PRECISION HIGH
C
C**********************************************************************
C     .. Scalar Arguments ..
      DOUBLE PRECISION HIGH,LOW
C     ..
C     .. External Functions ..
      DOUBLE PRECISION DRANFD
      EXTERNAL DRANFD
C     ..
C     .. Executable Statements ..
      IF (.NOT. (LOW.GT.HIGH)) GO TO 10
      WRITE (*,FMT=*) 'LOW > HIGH in DGENUNF: LOW ',LOW,' HIGH: ',HIGH
      WRITE (*,FMT=*) 'Abort'
      STOP 'LOW > High in DGENUNF - Abort'

   10 DGENUNF = LOW + (HIGH-LOW)*DRANFD()

      RETURN

      END
c -------------------------------------------------
      SUBROUTINE DGETCGN(G)
      INTEGER G
C**********************************************************************
C
C      SUBROUTINE DGETCGN(G)
C                         Get GeNerator
C
C     Returns in G the number of the current random number generator
C
C
C                              Arguments
C
C
C     G <-- Number of the current random number generator (1..32)
C                    INTEGER G
C
C**********************************************************************
C
      INTEGER CURNTG,NUMG
      SAVE CURNTG
      PARAMETER (NUMG=32)
      DATA CURNTG/1/
C
      G = CURNTG
      RETURN

      ENTRY SETCGN(G)
C**********************************************************************
C
C     SUBROUTINE SETCGN( G )
C                      Set GeNerator
C
C     Sets  the  current  generator to G.    All references to a
C     generat are to the current generator.
C
C
C                              Arguments
C
C
C     G --> Number of the current random number generator (1..32)
C                    INTEGER G
C
C**********************************************************************
C
C     Abort if generator number out of range
C
      IF (.NOT. (G.LT.0.OR.G.GT.NUMG)) GO TO 10
      WRITE (*,FMT=*) ' Generator number out of range in SETCGN:',
     +  ' Legal range is 1 to ',NUMG,' -- ABORT!'
      STOP ' Generator number out of range in SETCGN'

   10 CURNTG = G
      RETURN

      END
c -------------------------------------------------
      SUBROUTINE GETSD(ISEED1,ISEED2)
C**********************************************************************
C
C     SUBROUTINE GETSD(G,ISEED1,ISEED2)
C               GET SeeD
C
C     Returns the value of two integer seeds of the current generator
C
C     This is a transcription from Pascal to Fortran of routine
C     Get_State from the paper
C
C     L'Ecuyer, P. and  Cote,  S. "Implementing a Random Number
C     Package with Splitting Facilities."  ACM  Transactions on
C     Mathematical Software, 17:98-111 (1991)
C
C
C                              Arguments
C
C
C
C     ISEED1 <- First integer seed of generator G
C                                   INTEGER ISEED1
C
C     ISEED2 <- Second integer seed of generator G
C                                   INTEGER ISEED1
C
C**********************************************************************
C     .. Parameters ..
      INTEGER NUMG
      PARAMETER (NUMG=32)
C     ..
C     .. Scalar Arguments ..
      INTEGER ISEED1,ISEED2
C     ..
C     .. Scalars in Common ..
      INTEGER A1,A1VW,A1W,A2,A2VW,A2W,M1,M2
C     ..
C     .. Arrays in Common ..
      INTEGER CG1(NUMG),CG2(NUMG),IG1(NUMG),IG2(NUMG),LG1(NUMG),
     +        LG2(NUMG)
      LOGICAL QANTI(NUMG)
C     ..
C     .. Local Scalars ..
      INTEGER G
C     ..
C     .. External Functions ..
      LOGICAL QRGNIN
      EXTERNAL QRGNIN
C     ..
C     .. External Subroutines ..
      EXTERNAL DGETCGN
C     ..
C     .. Common blocks ..
      COMMON /GLOBE/M1,M2,A1,A2,A1W,A2W,A1VW,A2VW,IG1,IG2,LG1,LG2,CG1,
     +       CG2,QANTI
C     ..
C     .. Save statement ..
      SAVE /GLOBE/
C     ..
C     .. Executable Statements ..
C     Abort unless random number generator initialized
      IF (QRGNIN()) GO TO 10
      WRITE (*,FMT=*) ' GETSD called before random number generator ',
     +  ' initialized -- abort!'
      STOP ' GETSD called before random number generator initialized'

   10 CALL DGETCGN(G)
      ISEED1 = CG1(G)
      ISEED2 = CG2(G)
      RETURN

      END
c -------------------------------------------------
      INTEGER FUNCTION IGNLGI()
C**********************************************************************
C
C     INTEGER FUNCTION IGNLGI()
C               GeNerate LarGe Integer
C
C     Returns a random integer following a uniform distribution over
C     (1, 2147483562) using the current generator.
C
C     This is a transcription from Pascal to Fortran of routine
C     Random from the paper
C
C     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
C     with Splitting Facilities." ACM Transactions on Mathematical
C     Software, 17:98-111 (1991)
C
C**********************************************************************
C     .. Parameters ..
      INTEGER NUMG
      PARAMETER (NUMG=32)
C     ..
C     .. Scalars in Common ..
      INTEGER A1,A1VW,A1W,A2,A2VW,A2W,M1,M2
C     ..
C     .. Arrays in Common ..
      INTEGER CG1(NUMG),CG2(NUMG),IG1(NUMG),IG2(NUMG),LG1(NUMG),
     +        LG2(NUMG)
      LOGICAL QANTI(NUMG)
C     ..
C     .. Local Scalars ..
      INTEGER CURNTG,K,S1,S2,Z
      LOGICAL QQSSD
C     ..
C     .. External Functions ..
      LOGICAL QRGNIN
      EXTERNAL QRGNIN
C     ..
C     .. External Subroutines ..
      EXTERNAL DGETCGN,INRGCM,RGNQSD,SETALL
C     ..
C     .. Common blocks ..
      COMMON /GLOBE/M1,M2,A1,A2,A1W,A2W,A1VW,A2VW,IG1,IG2,LG1,LG2,CG1,
     +       CG2,QANTI
C     ..
C     .. Save statement ..
      SAVE /GLOBE/
C     ..
C     .. Executable Statements ..
C
C     IF THE RANDOM NUMBER PACKAGE HAS NOT BEEN INITIALIZED YET, DO SO.
C     IT CAN BE INITIALIZED IN ONE OF TWO WAYS : 1) THE FIRST CALL TO
C     THIS ROUTINE  2) A CALL TO SETALL.
C
      IF (.NOT. (QRGNIN())) CALL INRGCM
      CALL RGNQSD(QQSSD)
      IF (.NOT. (QQSSD)) CALL SETALL(1234567890,123456789)
C
C     Get Current Generator
C
      CALL DGETCGN(CURNTG)
      S1 = CG1(CURNTG)
      S2 = CG2(CURNTG)
      K = S1/53668
      S1 = A1* (S1-K*53668) - K*12211
      IF (S1.LT.0) S1 = S1 + M1
      K = S2/52774
      S2 = A2* (S2-K*52774) - K*3791
      IF (S2.LT.0) S2 = S2 + M2
      CG1(CURNTG) = S1
      CG2(CURNTG) = S2
      Z = S1 - S2
      IF (Z.LT.1) Z = Z + M1 - 1
      IF (QANTI(CURNTG)) Z = M1 - Z
      IGNLGI = Z
      RETURN

      END
c -------------------------------------------------
      INTEGER FUNCTION IGNUIN(LOW,HIGH)
C**********************************************************************
C
C     INTEGER FUNCTION IGNUIN( LOW, HIGH )
C
C               GeNerate Uniform INteger
C
C
C                              Function
C
C
C     Generates an integer uniformly distributed between LOW and HIGH.
C
C
C                              Arguments
C
C
C     LOW --> Low bound (inclusive) on integer value to be generated
C                         INTEGER LOW
C
C     HIGH --> High bound (inclusive) on integer value to be generated
C                         INTEGER HIGH
C
C
C                              Note
C
C
C     If (HIGH-LOW) > 2,147,483,561 prints error message on * unit and
C     stops the program.
C
C**********************************************************************

C     IGNLGI generates integers between 1 and 2147483562
C     MAXNUM is 1 less than maximum generable value
C     .. Parameters ..
      INTEGER MAXNUM
      PARAMETER (MAXNUM=2147483561)
      CHARACTER*(*) ERR1,ERR2
      PARAMETER (ERR1='LOW > HIGH in IGNUIN',
     +          ERR2=' ( HIGH - LOW ) > 2,147,483,561 in IGNUIN')
C     ..
C     .. Scalar Arguments ..
      INTEGER HIGH,LOW
C     ..
C     .. Local Scalars ..
      INTEGER ERR,IGN,MAXNOW,RANGE,RANP1
C     ..
C     .. External Functions ..
      INTEGER IGNLGI
      EXTERNAL IGNLGI
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC MOD
C     ..
C     .. Executable Statements ..
      IF (.NOT. (LOW.GT.HIGH)) GO TO 10
      ERR = 1
C      ABORT-PROGRAM
      GO TO 80

   10 RANGE = HIGH - LOW
      IF (.NOT. (RANGE.GT.MAXNUM)) GO TO 20
      ERR = 2
C      ABORT-PROGRAM
      GO TO 80

   20 IF (.NOT. (LOW.EQ.HIGH)) GO TO 30
      IGNUIN = LOW
      RETURN

      GO TO 70

C     Number to be generated should be in range 0..RANGE
C     Set MAXNOW so that the number of integers in 0..MAXNOW is an
C     integral multiple of the number in 0..RANGE

   30 RANP1 = RANGE + 1
      MAXNOW = (MAXNUM/RANP1)*RANP1
   40 IGN = IGNLGI() - 1
      IF (.NOT. (IGN.LE.MAXNOW)) GO TO 50
      IGNUIN = LOW + MOD(IGN,RANP1)
      RETURN

   50 GO TO 40

   60 CONTINUE
   70 CONTINUE
   80 IF (.NOT. (ERR.EQ.1)) GO TO 90
      WRITE (*,FMT=*) ERR1
      GO TO 100

C     TO ABORT-PROGRAM
   90 WRITE (*,FMT=*) ERR2
  100 WRITE (*,FMT=*) ' LOW: ',LOW,' HIGH: ',HIGH
      WRITE (*,FMT=*) ' Abort on Fatal ERROR'
      IF (.NOT. (ERR.EQ.1)) GO TO 110
      STOP 'LOW > HIGH in IGNUIN'

      GO TO 120

  110 STOP ' ( HIGH - LOW ) > 2,147,483,561 in IGNUIN'

  120 END
c -------------------------------------------------
      SUBROUTINE INITGN(ISDTYP)
C**********************************************************************
C
C     SUBROUTINE INITGN(ISDTYP)
C          INIT-ialize current G-e-N-erator
C
C     Reinitializes the state of the current generator
C
C     This is a transcription from Pascal to Fortran of routine
C     Init_Generator from the paper
C
C     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
C     with Splitting Facilities." ACM Transactions on Mathematical
C     Software, 17:98-111 (1991)
C
C
C                              Arguments
C
C
C     ISDTYP -> The state to which the generator is to be set
C
C          ISDTYP = -1  => sets the seeds to their initial value
C          ISDTYP =  0  => sets the seeds to the first value of
C                          the current block
C          ISDTYP =  1  => sets the seeds to the first value of
C                          the next block
C
C                                   INTEGER ISDTYP
C
C**********************************************************************
C     .. Parameters ..
      INTEGER NUMG
      PARAMETER (NUMG=32)
C     ..
C     .. Scalar Arguments ..
      INTEGER ISDTYP
C     ..
C     .. Scalars in Common ..
      INTEGER A1,A1VW,A1W,A2,A2VW,A2W,M1,M2
C     ..
C     .. Arrays in Common ..
      INTEGER CG1(NUMG),CG2(NUMG),IG1(NUMG),IG2(NUMG),LG1(NUMG),
     +        LG2(NUMG)
      LOGICAL QANTI(NUMG)
C     ..
C     .. Local Scalars ..
      INTEGER G
C     ..
C     .. External Functions ..
      LOGICAL QRGNIN
      INTEGER MLTMOD
      EXTERNAL QRGNIN,MLTMOD
C     ..
C     .. External Subroutines ..
      EXTERNAL DGETCGN
C     ..
C     .. Common blocks ..
      COMMON /GLOBE/M1,M2,A1,A2,A1W,A2W,A1VW,A2VW,IG1,IG2,LG1,LG2,CG1,
     +       CG2,QANTI
C     ..
C     .. Save statement ..
      SAVE /GLOBE/
C     ..
C     .. Executable Statements ..
C     Abort unless random number generator initialized
      IF (QRGNIN()) GO TO 10
      WRITE (*,FMT=*) ' INITGN called before random number generator ',
     +  ' initialized -- abort!'
      STOP ' INITGN called before random number generator initialized'

   10 CALL DGETCGN(G)
      IF ((-1).NE. (ISDTYP)) GO TO 20
      LG1(G) = IG1(G)
      LG2(G) = IG2(G)
      GO TO 50

   20 IF ((0).NE. (ISDTYP)) GO TO 30
      CONTINUE
      GO TO 50
C     do nothing
   30 IF ((1).NE. (ISDTYP)) GO TO 40
      LG1(G) = MLTMOD(A1W,LG1(G),M1)
      LG2(G) = MLTMOD(A2W,LG2(G),M2)
      GO TO 50

   40 STOP 'ISDTYP NOT IN RANGE'

   50 CG1(G) = LG1(G)
      CG2(G) = LG2(G)
      RETURN

      END
c -------------------------------------------------
      SUBROUTINE INRGCM
C**********************************************************************
C
C     SUBROUTINE INRGCM()
C          INitialize Random number Generator CoMmon
C
C
C                              Function
C
C
C     Initializes common area for random number generator. This saves
C     the nuisance of a BLOCK DATA routine and the difficulty of
C     assuring that the routine is loaded with the other routines.
C
C**********************************************************************
C     .. Parameters ..
      INTEGER NUMG
      PARAMETER (NUMG=32)
C     ..
C     .. Scalars in Common ..
      INTEGER A1,A1VW,A1W,A2,A2VW,A2W,M1,M2
C     ..
C     .. Arrays in Common ..
      INTEGER CG1(NUMG),CG2(NUMG),IG1(NUMG),IG2(NUMG),LG1(NUMG),
     +        LG2(NUMG)
      LOGICAL QANTI(NUMG)
C     ..
C     .. Local Scalars ..
      INTEGER I
      LOGICAL QDUM
C     ..
C     .. External Functions ..
      LOGICAL QRGNSN
      EXTERNAL QRGNSN
C     ..
C     .. Common blocks ..
      COMMON /GLOBE/M1,M2,A1,A2,A1W,A2W,A1VW,A2VW,IG1,IG2,LG1,LG2,CG1,
     +       CG2,QANTI
C     ..
C     .. Save statement ..
      SAVE /GLOBE/
C     ..
C     .. Executable Statements ..
C     V=20;                            W=30;
C
C     A1W = MOD(A1**(2**W),M1)         A2W = MOD(A2**(2**W),M2)
C     A1VW = MOD(A1**(2**(V+W)),M1)    A2VW = MOD(A2**(2**(V+W)),M2)
C
C   If V or W is changed A1W, A2W, A1VW, and A2VW need to be
C   recomputed.  An efficient way to precompute a**(2*j) MOD m is to
C   start with a and square it j times modulo m using the function
C   MLTMOD.
C
      M1 = 2147483563
      M2 = 2147483399
      A1 = 40014
      A2 = 40692
      A1W = 1033780774
      A2W = 1494757890
      A1VW = 2082007225
      A2VW = 784306273
      DO 10 I = 1,NUMG
          QANTI(I) = .FALSE.
   10 CONTINUE
C
C     Tell the world that common has been initialized
C
      QDUM = QRGNSN(.TRUE.)
      RETURN

      END
c -------------------------------------------------
      INTEGER FUNCTION MLTMOD(A,S,M)
C**********************************************************************
C
C     INTEGER FUNCTION MLTMOD(A,S,M)
C
C                    Returns (A*S) MOD M
C
C     This is a transcription from Pascal to Fortran of routine
C     MULtMod_Decompos from the paper
C
C     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
C     with Splitting Facilities." ACM Transactions on Mathematical
C     Software, 17:98-111 (1991)
C
C
C                              Arguments
C
C
C     A, S, M  -->
C                         INTEGER A,S,M
C
C**********************************************************************
C     .. Parameters ..
      INTEGER H
      PARAMETER (H=32768)
C     ..
C     .. Scalar Arguments ..
      INTEGER A,M,S
C     ..
C     .. Local Scalars ..
      INTEGER A0,A1,K,P,Q,QH,RH
C     ..
C     .. Executable Statements ..
C
C     H = 2**((b-2)/2) where b = 32 because we are using a 32 bit
C      machine. On a different machine recompute H
C
      IF (.NOT. (A.LE.0.OR.A.GE.M.OR.S.LE.0.OR.S.GE.M)) GO TO 10
      WRITE (*,FMT=*) ' A, M, S out of order in MLTMOD - ABORT!'
      WRITE (*,FMT=*) ' A = ',A,' S = ',S,' M = ',M
      WRITE (*,FMT=*) ' MLTMOD requires: 0 < A < M; 0 < S < M'
      STOP ' A, M, S out of order in MLTMOD - ABORT!'

   10 IF (.NOT. (A.LT.H)) GO TO 20
      A0 = A
      P = 0
      GO TO 120

   20 A1 = A/H
      A0 = A - H*A1
      QH = M/H
      RH = M - H*QH
      IF (.NOT. (A1.GE.H)) GO TO 50
      A1 = A1 - H
      K = S/QH
      P = H* (S-K*QH) - K*RH
   30 IF (.NOT. (P.LT.0)) GO TO 40
      P = P + M
      GO TO 30

   40 GO TO 60

   50 P = 0
C
C     P = (A2*S*H)MOD M
C
   60 IF (.NOT. (A1.NE.0)) GO TO 90
      Q = M/A1
      K = S/Q
      P = P - K* (M-A1*Q)
      IF (P.GT.0) P = P - M
      P = P + A1* (S-K*Q)
   70 IF (.NOT. (P.LT.0)) GO TO 80
      P = P + M
      GO TO 70

   80 CONTINUE
   90 K = P/QH
C
C     P = ((A2*H + A1)*S)MOD M
C
      P = H* (P-K*QH) - K*RH
  100 IF (.NOT. (P.LT.0)) GO TO 110
      P = P + M
      GO TO 100

  110 CONTINUE
  120 IF (.NOT. (A0.NE.0)) GO TO 150
C
C     P = ((A2*H + A1)*H*S)MOD M
C
      Q = M/A0
      K = S/Q
      P = P - K* (M-A0*Q)
      IF (P.GT.0) P = P - M
      P = P + A0* (S-K*Q)
  130 IF (.NOT. (P.LT.0)) GO TO 140
      P = P + M
      GO TO 130

  140 CONTINUE
  150 MLTMOD = P
C
      RETURN

      END
c -------------------------------------------------
      LOGICAL FUNCTION QRGNIN()
C**********************************************************************
C
C     LOGICAL FUNCTION QRGNIN()
C               Q Random GeNerators INitialized?
C
C     A trivial routine to determine whether or not the random
C     number generator has been initialized.  Returns .TRUE. if
C     it has, else .FALSE.
C
C**********************************************************************
C     .. Scalar Arguments ..
      LOGICAL QVALUE
C     ..
C     .. Local Scalars ..
      LOGICAL QINIT
C     ..
C     .. Entry Points ..
      LOGICAL QRGNSN
C     ..
C     .. Save statement ..
      SAVE QINIT
C     ..
C     .. Data statements ..
      DATA QINIT/.FALSE./
C     ..
C     .. Executable Statements ..
      QRGNIN = QINIT
      RETURN

      ENTRY QRGNSN(QVALUE)
C**********************************************************************
C
C     LOGICAL FUNCTION QRGNSN( QVALUE )
C               Q Random GeNerators Set whether iNitialized
C
C     Sets state of whether random number generator is initialized
C     to QVALUE.
C
C     This routine is actually an entry in QRGNIN, hence it is a
C     logical function.  It returns the (meaningless) value .TRUE.
C
C**********************************************************************
      QINIT = QVALUE
      QRGNSN = .TRUE.
      RETURN

      END
c -------------------------------------------------
      DOUBLE PRECISION FUNCTION DRANFD()
C**********************************************************************
C     NAME CHANGED BY DJS TO AVOID NAMING CONFLICTS
C**********************************************************************
C
C          DOUBLE PRECISION FUNCTION DRANFD()
C                RANDom number generator as a Function
C
C     Returns a random floating point number from a uniform
C     distribution over 0 - 1 (endpoints of this interval are not
C     returned) using the current generator
C
C     This is a transcription from Pascal to Fortran of routine
C     Uniform_01 from the paper
C
C     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
C     with Splitting Facilities." ACM Transactions on Mathematical
C     Software, 17:98-111 (1991)
C
C**********************************************************************
C     .. External Functions ..
      INTEGER IGNLGI
      EXTERNAL IGNLGI
C     ..
C     .. Executable Statements ..
C
C     4.656613057E-10 is 1/M1  M1 is set in a data statement in IGNLGI
C      and is currently 2147483563. If M1 changes, change this also.
C
      DRANFD = IGNLGI()*4.656613057D-10
      RETURN

      END
c -------------------------------------------------
      SUBROUTINE SETALL(ISEED1,ISEED2)
C**********************************************************************
C
C      SUBROUTINE SETALL(ISEED1,ISEED2)
C               SET ALL random number generators
C
C     Sets the initial seed of generator 1 to ISEED1 and ISEED2. The
C     initial seeds of the other generators are set accordingly, and
C     all generators states are set to these seeds.
C
C     This is a transcription from Pascal to Fortran of routine
C     Set_Initial_Seed from the paper
C
C     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
C     with Splitting Facilities." ACM Transactions on Mathematical
C     Software, 17:98-111 (1991)
C
C
C                              Arguments
C
C
C     ISEED1 -> First of two integer seeds
C                                   INTEGER ISEED1
C
C     ISEED2 -> Second of two integer seeds
C                                   INTEGER ISEED1
C
C**********************************************************************
C     .. Parameters ..
      INTEGER NUMG
      PARAMETER (NUMG=32)
C     ..
C     .. Scalar Arguments ..
      INTEGER ISEED1,ISEED2
      LOGICAL QSSD
C     ..
C     .. Scalars in Common ..
      INTEGER A1,A1VW,A1W,A2,A2VW,A2W,M1,M2
C     ..
C     .. Arrays in Common ..
      INTEGER CG1(NUMG),CG2(NUMG),IG1(NUMG),IG2(NUMG),LG1(NUMG),
     +        LG2(NUMG)
      LOGICAL QANTI(NUMG)
C     ..
C     .. Local Scalars ..
      INTEGER G,OCGN
      LOGICAL QQSSD
C     ..
C     .. External Functions ..
      INTEGER MLTMOD
      LOGICAL QRGNIN
      EXTERNAL MLTMOD,QRGNIN
C     ..
C     .. External Subroutines ..
      EXTERNAL DGETCGN,INITGN,INRGCM,SETCGN
C     ..
C     .. Common blocks ..
      COMMON /GLOBE/M1,M2,A1,A2,A1W,A2W,A1VW,A2VW,IG1,IG2,LG1,LG2,CG1,
     +       CG2,QANTI
C     ..
C     .. Save statement ..
      SAVE /GLOBE/,QQSSD
C     ..
C     .. Data statements ..
      DATA QQSSD/.FALSE./
C     ..
C     .. Executable Statements ..
C
C     TELL IGNLGI, THE ACTUAL NUMBER GENERATOR, THAT THIS ROUTINE
C      HAS BEEN CALLED.
C
      QQSSD = .TRUE.
      CALL DGETCGN(OCGN)
C
C     Initialize Common Block if Necessary
C
      IF (.NOT. (QRGNIN())) CALL INRGCM
      IG1(1) = ISEED1
      IG2(1) = ISEED2
      CALL INITGN(-1)
      DO 10 G = 2,NUMG
          IG1(G) = MLTMOD(A1VW,IG1(G-1),M1)
          IG2(G) = MLTMOD(A2VW,IG2(G-1),M2)
          CALL SETCGN(G)
          CALL INITGN(-1)
   10 CONTINUE
      CALL SETCGN(OCGN)
      RETURN

      ENTRY RGNQSD(QSSD)
C**********************************************************************
C
C     SUBROUTINE RGNQSD
C                    Random Number Generator Query SeeD set?
C
C     Returns (LOGICAL) QSSD as .TRUE. if SETALL has been invoked,
C     otherwise returns .FALSE.
C
C**********************************************************************
      QSSD = QQSSD
      RETURN

      END
c -------------------------------------------------
      SUBROUTINE SETSD(ISEED1,ISEED2)
C**********************************************************************
C
C     SUBROUTINE SETSD(ISEED1,ISEED2)
C               SET S-ee-D of current generator
C
C     Resets the initial seed of the current generator to ISEED1 and
C     ISEED2. The seeds of the other generators remain unchanged.
C
C     This is a transcription from Pascal to Fortran of routine
C     Set_Seed from the paper
C
C     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
C     with Splitting Facilities." ACM Transactions on Mathematical
C     Software, 17:98-111 (1991)
C
C
C                              Arguments
C
C
C     ISEED1 -> First integer seed
C                                   INTEGER ISEED1
C
C     ISEED2 -> Second integer seed
C                                   INTEGER ISEED1
C
C**********************************************************************
C     .. Parameters ..
      INTEGER NUMG
      PARAMETER (NUMG=32)
C     ..
C     .. Scalar Arguments ..
      INTEGER ISEED1,ISEED2
C     ..
C     .. Scalars in Common ..
      INTEGER A1,A1VW,A1W,A2,A2VW,A2W,M1,M2
C     ..
C     .. Arrays in Common ..
      INTEGER CG1(NUMG),CG2(NUMG),IG1(NUMG),IG2(NUMG),LG1(NUMG),
     +        LG2(NUMG)
      LOGICAL QANTI(NUMG)
C     ..
C     .. Local Scalars ..
      INTEGER G
C     ..
C     .. External Functions ..
      LOGICAL QRGNIN
      EXTERNAL QRGNIN
C     ..
C     .. External Subroutines ..
      EXTERNAL DGETCGN,INITGN
C     ..
C     .. Common blocks ..
      COMMON /GLOBE/M1,M2,A1,A2,A1W,A2W,A1VW,A2VW,IG1,IG2,LG1,LG2,CG1,
     +       CG2,QANTI
C     ..
C     .. Save statement ..
      SAVE /GLOBE/
C     ..
C     .. Executable Statements ..
C     Abort unless random number generator initialized
      IF (QRGNIN()) GO TO 10
      WRITE (*,FMT=*) ' SETSD called before random number generator ',
     +  ' initialized -- abort!'
      STOP ' SETSD called before random number generator initialized'

   10 CALL DGETCGN(G)
      IG1(G) = ISEED1
      IG2(G) = ISEED2
      CALL INITGN(-1)
      RETURN

      END
c -------------------------------------------------
      DOUBLE PRECISION FUNCTION DSEXPO()
      DOUBLE PRECISION Q
      DOUBLE PRECISION Q1
      DOUBLE PRECISION A
      DOUBLE PRECISION U
      DOUBLE PRECISION DRANFD
      DOUBLE PRECISION USTAR
      DOUBLE PRECISION UMIN
C*********************************************************************C
C                                                                     C
C                                                                     C
C    (STANDARD-)  E X P O N E N T I A L   DISTRIBUTION                C
C                                                                     C
C                                                                     C
C*********************************************************************C
**********************************************************************C
C                                                                     C
C    FOR DETAILS SEE:                                                 C
C                                                                     C
C              AHRENS, J.H. AND DIETER, U.                            C
C              COMPUTER METHODS FOR SAMPLING FROM THE                 C
C              EXPONENTIAL AND NORMAL DISTRIBUTIONS.                  C
C              COMM. ACM, 15,10 (OCT. 1972), 873 - 882.               C
C                                                                     C
C    ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM       C
C    'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)       C
C                                                                     C
C    Modified by Barry W. Brown, Feb 3, 1988 to use DRANFD instead of C
C    SUNIF.  The argument IR thus goes away.                          C
C                                                                     C
C*********************************************************************C
C
      DIMENSION Q(8)
      EQUIVALENCE (Q(1),Q1)
C
C     Q(N) = SUM(ALOG(2.0)**K/K!)    K=1,..,N ,      THE HIGHEST N
C     (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
C
      DATA Q/.6931472D0,.9333737D0,.9888778D0,.9984959D0,.9998293D0,
     +     .9999833D0,.9999986D0,.9999999D0/
C
   10 A = 0.0D0
      U = DRANFD()
      GO TO 30

   20 A = A + Q1
   30 U = U + U
      IF (U.LE.1.0D0) GO TO 20
   40 U = U - 1.0D0
      IF (U.GT.Q1) GO TO 60
   50 DSEXPO = A + U
      RETURN

   60 I = 1
      USTAR = DRANFD()
      UMIN = USTAR
   70 USTAR = DRANFD()
      IF (USTAR.LT.UMIN) UMIN = USTAR
   80 I = I + 1
      IF (U.GT.Q(I)) GO TO 70
   90 DSEXPO = A + UMIN*Q1
      RETURN

      END
c -------------------------------------------------
      DOUBLE PRECISION FUNCTION DSGAMMA(A)
      DOUBLE PRECISION A
      DOUBLE PRECISION Q1
      DOUBLE PRECISION Q2
      DOUBLE PRECISION Q3
      DOUBLE PRECISION Q4
      DOUBLE PRECISION Q5
      DOUBLE PRECISION Q6
      DOUBLE PRECISION Q7
      DOUBLE PRECISION A1
      DOUBLE PRECISION A2
      DOUBLE PRECISION A3
      DOUBLE PRECISION A4
      DOUBLE PRECISION A5
      DOUBLE PRECISION A6
      DOUBLE PRECISION A7
      DOUBLE PRECISION E1
      DOUBLE PRECISION E2
      DOUBLE PRECISION E3
      DOUBLE PRECISION E4
      DOUBLE PRECISION E5
      DOUBLE PRECISION AA
      DOUBLE PRECISION AAA
      DOUBLE PRECISION SQRT32
      DOUBLE PRECISION S2
      DOUBLE PRECISION S
      DOUBLE PRECISION D
      DOUBLE PRECISION Q0
      DOUBLE PRECISION B
      DOUBLE PRECISION SI
      DOUBLE PRECISION C
      DOUBLE PRECISION T
      DOUBLE PRECISION DSNORM
      DOUBLE PRECISION X
      DOUBLE PRECISION U
      DOUBLE PRECISION DRANFD
      DOUBLE PRECISION R
      DOUBLE PRECISION V
      DOUBLE PRECISION Q
      DOUBLE PRECISION E
      DOUBLE PRECISION DSEXPO
      DOUBLE PRECISION W
      DOUBLE PRECISION P
C*********************************************************************C
C                                                                     C
C                                                                     C
C    (STANDARD-)  G A M M A  DISTRIBUTION                             C
C                                                                     C
C                                                                     C
C*********************************************************************C
C*********************************************************************C
C                                                                     C
C              PARAMETER  A >= 1.0  !                                 C
C                                                                     C
C*********************************************************************C
C                                                                     C
C    FOR DETAILS SEE:                                                 C
C                                                                     C
C              AHRENS, J.H. AND DIETER, U.                            C
C              GENERATING GAMMA VARIATES BY A                         C
C              MODIFIED REJECTION TECHNIQUE.                          C
C              COMM. ACM, 25,1 (JAN. 1982), 47 - 54.                  C
C                                                                     C
C    STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER     C
C                                (STRAIGHTFORWARD IMPLEMENTATION)     C
C                                                                     C
C    Modified by Barry W. Brown, Feb 3, 1988 to use DRANFD instead of C
C    SUNIF.  The argument IR thus goes away.                          C
C                                                                     C
C*********************************************************************C
C                                                                     C
C              PARAMETER  0.0 < A < 1.0  !                            C
C                                                                     C
C*********************************************************************C
C                                                                     C
C    FOR DETAILS SEE:                                                 C
C                                                                     C
C              AHRENS, J.H. AND DIETER, U.                            C
C              COMPUTER METHODS FOR SAMPLING FROM GAMMA,              C
C              BETA, POISSON AND BINOMIAL DISTRIBUTIONS.              C
C              COMPUTING, 12 (1974), 223 - 246.                       C
C                                                                     C
C    (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER)    C
C                                                                     C
C*********************************************************************C
C
C    INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION
C    OUTPUT: DSGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION
C
C    COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K))
C    COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K)
C    COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K)
C
      DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7/.04166669D0,.02083148D0,.00801191D0,
     +     .00144121D0,-.00007388D0,.00024511D0,.00024240D0/
      DATA A1,A2,A3,A4,A5,A6,A7/.3333333D0,-.2500030D0,.2000062D0,
     +     -.1662921D0,.1423657D0,-.1367177D0,.1233795D0/
      DATA E1,E2,E3,E4,E5/1.D0,.4999897D0,.1668290D0,.0407753D0,
     +     .0102930D0/
C
C     PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A"
C     SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380
C
      DATA AA/0.0D0/,AAA/0.0D0/,SQRT32/5.656854D0/
C
C     SAVE STATEMENTS
C
      SAVE AA,AAA,S2,S,D,Q0,B,SI,C
C
      IF (A.EQ.AA) GO TO 10
      IF (A.LT.1.0D0) GO TO 140
C
C     STEP  1:  RECALCULATIONS OF S2,S,D IF A HAS CHANGED
C
      AA = A
      S2 = A - 0.5D0
      S = SQRT(S2)
      D = SQRT32 - 12.0D0*S
C
C     STEP  2:  T=STANDARD NORMAL DEVIATE,
C               X=(S,1/2)-NORMAL DEVIATE.
C               IMMEDIATE ACCEPTANCE (I)
C
   10 T = DSNORM()
      X = S + 0.5D0*T
      DSGAMMA = X*X
      IF (T.GE.0.0D0) RETURN
C
C     STEP  3:  U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S)
C
      U = DRANFD()
      IF (D*U.LE.T*T*T) RETURN
C
C     STEP  4:  RECALCULATIONS OF Q0,B,SI,C IF NECESSARY
C
      IF (A.EQ.AAA) GO TO 40
      AAA = A
      R = 1.0D0/A
      Q0 = ((((((Q7*R+Q6)*R+Q5)*R+Q4)*R+Q3)*R+Q2)*R+Q1)*R
C
C               APPROXIMATION DEPENDING ON SIZE OF PARAMETER A
C               THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND
C               C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS
C
      IF (A.LE.3.686D0) GO TO 30
      IF (A.LE.13.022D0) GO TO 20
C
C               CASE 3:  A .GT. 13.022
C
      B = 1.77D0
      SI = .75D0
      C = .1515D0/S
      GO TO 40
C
C               CASE 2:  3.686 .LT. A .LE. 13.022
C
   20 B = 1.654D0 + .0076D0*S2
      SI = 1.68D0/S + .275D0
      C = .062D0/S + .024D0
      GO TO 40
C
C               CASE 1:  A .LE. 3.686
C
   30 B = .463D0 + S + .178D0*S2
      SI = 1.235D0
      C = .195D0/S - .079D0 + .16D0*S
C
C     STEP  5:  NO QUOTIENT TEST IF X NOT POSITIVE
C
   40 IF (X.LE.0.0D0) GO TO 70
C
C     STEP  6:  CALCULATION OF V AND QUOTIENT Q
C
      V = T/ (S+S)
      IF (ABS(V).LE.0.25D0) GO TO 50
      Q = Q0 - S*T + 0.25D0*T*T + (S2+S2)*DLOG(1.0D0+V)
      GO TO 60

   50 Q=Q0+0.5D0*T*T*((((((A7*V+A6)*V+A5)*V+A4)*V+A3)*V+A2)*V+A1)*V
C
C     STEP  7:  QUOTIENT ACCEPTANCE (Q)
C
   60 IF (DLOG(1.0D0-U).LE.Q) RETURN
C
C     STEP  8:  E=STANDARD EXPONENTIAL DEVIATE
C               U= 0,1 -UNIFORM DEVIATE
C               T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE
C
   70 E = DSEXPO()
      U = DRANFD()
      U = U + U - 1.0D0
      T = B + SIGN(SI*E,U)
      IF (.NOT. (U.GE.0.0D0)) GO TO 80
      T = B + SI*E
      GO TO 90

   80 T = B - SI*E

C
C     STEP  9:  REJECTION IF T .LT. TAU(1) = -.71874483771719
C
   90 IF (T.LT. (-.7187449D0)) GO TO 70
C
C     STEP 10:  CALCULATION OF V AND QUOTIENT Q
C
      V = T/ (S+S)
      IF (ABS(V).LE.0.25D0) GO TO 100
      Q = Q0 - S*T + 0.25D0*T*T + (S2+S2)*DLOG(1.0D0+V)
      GO TO 110

  100 Q=Q0+0.5D0*T*T*((((((A7*V+A6)*V+A5)*V+A4)*V+A3)*V+A2)*V+A1)*V
C
C     STEP 11:  HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8)
C
  110 IF (Q.LE.0.0D0) GO TO 70
      IF (Q.LE.0.5D0) GO TO 120
      W = EXP(Q) - 1.0D0
      GO TO 130

  120 W = ((((E5*Q+E4)*Q+E3)*Q+E2)*Q+E1)*Q
C
C               IF T IS REJECTED, SAMPLE AGAIN AT STEP 8
C
  130 IF (C*ABS(U).GT.W*EXP(E-0.5D0*T*T)) GO TO 70
      X = S + 0.5D0*T
      DSGAMMA = X*X
      RETURN
C
C     ALTERNATE METHOD FOR PARAMETERS A BELOW 1  (.3678794=EXP(-1.))
C
  140 AA = 0.0D0
      B = 1.0D0 + .3678794D0*A
  150 P = B*DRANFD()
      IF (P.GE.1.0D0) GO TO 160
      DSGAMMA = EXP(DLOG(P)/A)
      IF (DSEXPO().LT.DSGAMMA) GO TO 150
      RETURN

  160 DSGAMMA = -DLOG((B-P)/A)
      IF (DSEXPO().LT. (1.0D0-A)*DLOG(DSGAMMA)) GO TO 150
      RETURN

      END
c -------------------------------------------------
      DOUBLE PRECISION FUNCTION DSNORM()
      DOUBLE PRECISION A
      DOUBLE PRECISION D
      DOUBLE PRECISION T
      DOUBLE PRECISION H
      DOUBLE PRECISION U
      DOUBLE PRECISION DRANFD
      DOUBLE PRECISION S
      DOUBLE PRECISION USTAR
      DOUBLE PRECISION AA
      DOUBLE PRECISION W
      DOUBLE PRECISION Y
      DOUBLE PRECISION TT
C*********************************************************************C
C                                                                     C
C                                                                     C
C    (STANDARD-)  N O R M A L  DISTRIBUTION                           C
C                                                                     C
C                                                                     C
C*********************************************************************C
C*********************************************************************C
C                                                                     C
C    FOR DETAILS SEE:                                                 C
C                                                                     C
C              AHRENS, J.H. AND DIETER, U.                            C
C              EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM             C
C              SAMPLING FROM THE NORMAL DISTRIBUTION.                 C
C              MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937.          C
C                                                                     C
C    ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL'  C
C    (M=5) IN THE ABOVE PAPER     (SLIGHTLY MODIFIED IMPLEMENTATION)  C
C                                                                     C
C    Modified by Barry W. Brown, Feb 3, 1988 to use DRANFD instead of C
C    SUNIF.  The argument IR thus goes away.                          C
C                                                                     C
C*********************************************************************C
C
      DIMENSION A(32),D(31),T(31),H(31)
C
C     THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
C     H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
C
      DATA A/0.0D0,.3917609D-1,.7841241D-1,.1177699D0,.1573107D0,
     +     .1970991D0,.2372021D0,.2776904D0,.3186394D0,.3601299D0,
     +     .4022501D0,.4450965D0,.4887764D0,.5334097D0,.5791322D0,
     +     .6260990D0,.6744898D0,.7245144D0,.7764218D0,.8305109D0,
     +     .8871466D0,.9467818D0,1.009990D0,1.077516D0,1.150349D0,
     +     1.229859D0,1.318011D0,1.417797D0,1.534121D0,1.675940D0,
     +     1.862732D0,2.153875D0/
      DATA D/5*0.0D0,.2636843D0,.2425085D0,.2255674D0,.2116342D0,
     +     .1999243D0,.1899108D0,.1812252D0,.1736014D0,.1668419D0,
     +     .1607967D0,.1553497D0,.1504094D0,.1459026D0,.1417700D0,
     +     .1379632D0,.1344418D0,.1311722D0,.1281260D0,.1252791D0,
     +     .1226109D0,.1201036D0,.1177417D0,.1155119D0,.1134023D0,
     +     .1114027D0,.1095039D0/
      DATA T/.7673828D-3,.2306870D-2,.3860618D-2,.5438454D-2,
     +     .7050699D-2,.8708396D-2,.1042357D-1,.1220953D-1,.1408125D-1,
     +     .1605579D-1,.1815290D-1,.2039573D-1,.2281177D-1,.2543407D-1,
     +     .2830296D-1,.3146822D-1,.3499233D-1,.3895483D-1,.4345878D-1,
     +     .4864035D-1,.5468334D-1,.6184222D-1,.7047983D-1,.8113195D-1,
     +     .9462444D-1,.1123001D0,.1364980D0,.1716886D0,.2276241D0,
     +     .3304980D0,.5847031D0/
      DATA H/.3920617D-1,.3932705D-1,.3950999D-1,.3975703D-1,
     +     .4007093D-1,.4045533D-1,.4091481D-1,.4145507D-1,.4208311D-1,
     +     .4280748D-1,.4363863D-1,.4458932D-1,.4567523D-1,.4691571D-1,
     +     .4833487D-1,.4996298D-1,.5183859D-1,.5401138D-1,.5654656D-1,
     +     .5953130D-1,.6308489D-1,.6737503D-1,.7264544D-1,.7926471D-1,
     +     .8781922D-1,.9930398D-1,.1155599D0,.1404344D0,.1836142D0,
     +     .2790016D0,.7010474D0/
C
   10 U = DRANFD()
      S = 0.0D0
      IF (U.GT.0.5D0) S = 1.0D0
      U = U + U - S
   20 U = 32.0D0*U
      I = INT(U)
      IF (I.EQ.32) I = 31
      IF (I.EQ.0) GO TO 100
C
C                                START CENTER
C
   30 USTAR = U - DBLE(I)
      AA = A(I)
   40 IF (USTAR.LE.T(I)) GO TO 60
      W = (USTAR-T(I))*H(I)
C
C                                EXIT   (BOTH CASES)
C
   50 Y = AA + W
      DSNORM = Y
      IF (S.EQ.1.0D0) DSNORM = -Y
      RETURN
C
C                                CENTER CONTINUED
C
   60 U = DRANFD()
      W = U* (A(I+1)-AA)
      TT = (0.5D0*W+AA)*W
      GO TO 80

   70 TT = U
      USTAR = DRANFD()
   80 IF (USTAR.GT.TT) GO TO 50
   90 U = DRANFD()
      IF (USTAR.GE.U) GO TO 70
      USTAR = DRANFD()
      GO TO 40
C
C                                START TAIL
C
  100 I = 6
      AA = A(32)
      GO TO 120

  110 AA = AA + D(I)
      I = I + 1
  120 U = U + U
      IF (U.LT.1.0D0) GO TO 110
  130 U = U - 1.0D0
  140 W = U*D(I)
      TT = (0.5D0*W+AA)*W
      GO TO 160

  150 TT = U
  160 USTAR = DRANFD()
      IF (USTAR.GT.TT) GO TO 50
  170 U = DRANFD()
      IF (USTAR.GE.U) GO TO 150
      U = DRANFD()
      GO TO 140

      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
      SUBROUTINE ABORT                                                  MLu
C
C     TEMPORARY ABORT ROUTINE FOR NON-HARRIS SITES
C
      WRITE (6,10)                                                      MLu
 10   FORMAT (//,' ***** PROGRAM CALLED ABORT *****',//)                MLu
C
C     CALL EXIT (2)                                                     u
C
C     GET THE SQUARE ROOT OF A NEGATIVE NUMBER TO CAUSE A FATAL ERROR
      X = SIN (0.2)                                                     L
      X = -X                                                            L
      Y = SQRT (X)                                                      L
      WRITE (*,*) Y                                                     L
      STOP                                                              MLu
      END                                                               MLu

max_stars_repo_name: HydrologicEngineeringCenter/heclib

max_issues_repo_name: HydrologicEngineeringCenter/heclib

max_forks_repo_name: HydrologicEngineeringCenter/heclib


text:
program ramses2tipsy

  use random
  use io_ramses

  implicit none

  integer::ndim,n,i,j,k,iii,twotondim,ncoarse,indcell
  integer::ivar,nvar,ncpu,lmax=100,levelmin,iskip
  integer::nx,ny,nz
  integer::nlevelmax,ilevel1,ngrid1
  integer::nlevelmaxs,nlevel,iout
  integer::ipos
  integer::ngridmax,nstep_coarse,icpu
  integer::nhx,nhy,ihx,ihy,ivar1,ivar2
  real(KIND=8)::gamma,smallr,smallc,gammah
  real(KIND=8)::boxlen,boxlen2
  real(KIND=8)::t,aexp,h0,t2,aexp2,hexp2
  real(KIND=8)::omega_m,omega_l,omega_k,omega_b
  real(kind=8)::omega_m2,omega_l2,omega_k2,omega_b2
  real(kind=8)::scale_l,scale_d,scale_t
  real(kind=8)::metmax=0d0

  integer::nx_sample=0,ny_sample=0,ngridtot
  integer::ngrid,imin,imax,jmin,jmax,kmin,kmax
  integer::ncpu2,npart2,ndim2,nlevelmax2,nstep_coarse2
  integer::nx2,ny2,nz2,ngridmax2,nvarh,ndimh,nlevelmaxh
  integer::nx_full,ny_full,lmin,nboundary,ngrid_current,lllmin
  integer::ix,iy,iz,ndom,impi,bit_length,maxdom,ii,jj,kk
  integer,dimension(1:8)::idom,jdom,kdom,cpu_min,cpu_max
  real(KIND=8),dimension(1:8)::bounding_min,bounding_max
  real(KIND=8)::dkey,order_min,dmax,ddx,dxline,ddy,dex,dey,weight,msph=0.d0
  real(KIND=8)::xmin=0,xmax=-1,ymin=0,ymax=-1,zmin=0,zmax=-1,mdm=0.d0,mres=0.d0

  integer::ilevel,ncpu_read,three
  real(kind=8)deltax
  character(len=3)::typ='all'
  character(LEN=5)::nchar
  character(LEN=80)::ordering
  character(LEN=128)::nomfich
  character(LEN=128)::repository,outfich,filetype='bin'
  character(LEN=13)::string
  logical::ok,ok_part,ok_cell,do_max,do_id=.false.
  real(kind=8),dimension(:),allocatable::bound_key,xdp
  logical,dimension(:),allocatable::cpu_read
  integer,dimension(:),allocatable::cpu_list

  integer::ndummypart,nmin=0,nmax=0,nold,nnold,ndummyold
  integer::partcount,respart,denspartcount
  real(KIND=8)::dummy,partmass,volume,facdens=0.d0,averdens
  integer::delm=0,levelsel

  integer::npart,nstar_tot,nsink,npart_tot,nsink_tot
  integer::npart_actual,ndm_actual,nstar_actual,npartnow
  character(LEN=5)::nn
  logical::hydrok=.false.,partok=.false.,metal=.false.,metgas=.false.
  logical::cosmo=.false.,star=.false.,sink=.false.,mhd=.false.,gas=.true.

  integer::n_frw
  real(KIND=8)::time,time_tot,time_simu,time_uni
  real(KIND=8),dimension(:),allocatable::aexp_frw,hexp_frw,tau_frw,t_frw

  integer ,dimension(1:1,1:IRandNumSize)::allseed
  integer ,dimension(1:IRandNumSize)::localseed
  integer::iseed=0,poisson

  call read_params

  !-----------------------------------------------
  ! Reading files in RAMSES format
  !-----------------------------------------------
  ipos=INDEX(repository,'output_')
  nchar=repository(ipos+7:ipos+13)
  nomfich=TRIM(repository)//'/hydro_'//TRIM(nchar)//'.out00001'
  inquire(file=nomfich, exist=ok) ! verify input file
  hydrok=ok
  if ( .not. ok ) then
     print *,TRIM(nomfich)//' not found.'
  endif
  nomfich=TRIM(repository)//'/part_'//TRIM(nchar)//'.out00001'
  inquire(file=nomfich, exist=ok) ! verify input file
  partok=ok
  if ( .not. ok ) then
     print *,TRIM(nomfich)//' not found.'
  endif
  nomfich=TRIM(repository)//'/amr_'//TRIM(nchar)//'.out00001'
  inquire(file=nomfich, exist=ok) ! verify input file
  if ( .not. ok ) then
     print *,TRIM(nomfich)//' not found.'
     stop
  endif

  nomfich=TRIM(repository)//'/amr_'//TRIM(nchar)//'.out00001'
  open(unit=10,file=nomfich,status='old',form='unformatted')
  read(10)ncpu
  read(10)ndim
  read(10)nx,ny,nz
  read(10)nlevelmax
  read(10)ngridmax
  read(10)nboundary
  read(10)ngrid_current
  read(10)boxlen
  read(10)
  read(10)
  read(10)
  read(10)
  read(10)
  read(10)
  read(10)
  read(10)
  read(10)
  read(10)
  read(10)msph
  close(10)
  twotondim=2**ndim

  ! Rescaling mass_sph
  msph=msph/boxlen**3

  ! Default values for box size
  if(xmax<0)xmax=boxlen
  if(ymax<0)ymax=boxlen
  if(zmax<0)zmax=boxlen

  if(hydrok)then
     ! Read nvarh from the Hydro file
     nomfich=TRIM(repository)//'/hydro_'//TRIM(nchar)//'.out00001'
     open(unit=12,file=nomfich,status='old',form='unformatted')
     read(12)
     read(12)nvarh
     read(12)
     read(12)
     read(12)
     read(12)gamma
     close(12)
     if(ndim==2)then
        if((.not.mhd).and.(nvarh<=4))then
           metgas=.false.
           metal=.false.
        end if
        if((.not.mhd).and.(nvarh>4))then
           metgas=.true.
        end if
        if((mhd).and.(nvarh<=11))then
           metgas=.false.
           metal=.false.
        end if
        if((mhd).and.(nvarh>11))then
           metgas=.true.
        end if
     end if
     if(ndim==3)then
        if((.not.mhd).and.(nvarh<=5))then
           metgas=.false.
           metal=.false.
        end if
        if((.not.mhd).and.(nvarh>5))then
           metgas=.true.
        end if
        if((mhd).and.(nvarh<=11))then
           metgas=.false.
           metal=.false.
        end if
        if((mhd).and.(nvarh>11))then
           metgas=.true.
        end if
     end if
  end if

  if(partok)then
     npart_tot=0
     nsink_tot=0
     do i=1,ncpu
        ! Read number of particles from the Part file
        write(nn,'(I5.5)')i
        nomfich=TRIM(repository)//'/part_'//TRIM(nchar)//'.out'//nn
        open(unit=11,file=nomfich,status='old',form='unformatted')
        read(11)
        read(11)
        read(11)npart
        read(11)
        read(11)nstar_tot
        read(11)
        read(11)
        read(11)nsink
        close(11)
        nsink_tot=nsink_tot+nsink
        npart_tot=npart_tot+npart
     enddo
     if(nsink_tot>0)sink=.true.
     if(nstar_tot>0)star=.true.
     if(metgas.and.star)metal=.true.
  else
     sink=.false.
     star=.false.
     metal=.false.
  end if

  nomfich=TRIM(repository)//'/info_'//TRIM(nchar)//'.txt'
  inquire(file=nomfich, exist=ok) ! verify input file
  if ( .not. ok ) then
     print *,TRIM(nomfich)//' not found.'
     stop
  endif
  open(unit=10,file=nomfich,form='formatted',status='old')
  read(10,*)
  read(10,*)
  read(10,'(A13,I11)')string,levelmin
  read(10,*)
  read(10,*)
  read(10,*)
  read(10,*)
  read(10,*)
  read(10,'(A13,E23.15)')string,t
  read(10,'(A13,E23.15)')string,aexp
  read(10,'(A13,E23.15)')string,h0
  read(10,'(A13,E23.15)')string,omega_m
  read(10,'(A13,E23.15)')string,omega_l
  read(10,'(A13,E23.15)')string,omega_k
  read(10,'(A13,E23.15)')string,omega_b
  read(10,'(A13,E23.15)')string,scale_l
  read(10,'(A13,E23.15)')string,scale_d
  read(10,'(A13,E23.15)')string,scale_t
  read(10,*)
  read(10,'("ordering type=",A80)'),ordering
  read(10,*)
  allocate(cpu_list(1:ncpu))
  if(TRIM(ordering).eq.'hilbert')then
     allocate(bound_key(0:ncpu))
     allocate(cpu_read(1:ncpu))
     cpu_read=.false.
     do impi=1,ncpu
        read(10,'(I8,1X,E23.15,1X,E23.15)')i,bound_key(impi-1),bound_key(impi)
     end do
  endif
  close(10)

  lmax=max(min(lmax,nlevelmax),1)
  levelsel=levelmin+delm

  if(TRIM(ordering).eq.'hilbert')then

     dmax=max(xmax-xmin,ymax-ymin,zmax-zmin)
     do ilevel=1,lmax
        deltax=0.5d0**ilevel
        if(deltax.lt.dmax)exit
     end do
     lllmin=ilevel
     bit_length=lllmin-1
     maxdom=2**bit_length
     imin=0; imax=0; jmin=0; jmax=0; kmin=0; kmax=0
     if(bit_length>0)then
        imin=int(xmin*dble(maxdom))
        imax=imin+1
        jmin=int(ymin*dble(maxdom))
        jmax=jmin+1
        kmin=int(zmin*dble(maxdom))
        kmax=kmin+1
     endif

     dkey=(dble(2**(lmax+1)/dble(maxdom)))**ndim

     ndom=1
     if(bit_length>0)ndom=8
     idom(1)=imin; idom(2)=imax
     idom(3)=imin; idom(4)=imax
     idom(5)=imin; idom(6)=imax
     idom(7)=imin; idom(8)=imax
     jdom(1)=jmin; jdom(2)=jmin
     jdom(3)=jmax; jdom(4)=jmax
     jdom(5)=jmin; jdom(6)=jmin
     jdom(7)=jmax; jdom(8)=jmax
     kdom(1)=kmin; kdom(2)=kmin
     kdom(3)=kmin; kdom(4)=kmin
     kdom(5)=kmax; kdom(6)=kmax
     kdom(7)=kmax; kdom(8)=kmax

     do i=1,ndom
        if(bit_length>0)then
           ii=idom(i)
           jj=jdom(i)
           kk=kdom(i)
           call hilbert3d(ii,jj,kk,order_min,bit_length,1)
        else
           order_min=0.0d0
        endif
        bounding_min(i)=(order_min)*dkey
        bounding_max(i)=(order_min+1.0d0)*dkey
     end do
     cpu_min=0; cpu_max=0
     do impi=1,ncpu
        do i=1,ndom
           if (bound_key(impi-1).le.bounding_min(i).and.&
                &bound_key(impi).gt.bounding_min(i))then
              cpu_min(i)=impi
           endif
           if (bound_key(impi-1).lt.bounding_max(i).and.&
                &bound_key(impi).ge.bounding_max(i))then
              cpu_max(i)=impi
           endif
        end do
     end do

     ncpu_read=0
     do i=1,ndom
        do j=cpu_min(i),cpu_max(i)
           if(.not. cpu_read(j))then
              ncpu_read=ncpu_read+1
              cpu_list(ncpu_read)=j
              cpu_read(j)=.true.
           endif
        enddo
     enddo
  else
     ncpu_read=ncpu
     do j=1,ncpu
        cpu_list(j)=j
     end do
  end  if

  if(cosmo)then
     !-----------------------
     ! Cosmological model
     !-----------------------
     ! Allocate look-up tables
     n_frw=1000
     allocate(aexp_frw(0:n_frw),hexp_frw(0:n_frw))
     allocate(tau_frw(0:n_frw),t_frw(0:n_frw))

     ! Compute Friedman model look up table
     write(*,*)'Computing Friedman model'
     call friedman(dble(omega_m),dble(omega_l),dble(omega_k), &
          & 1.d-6,1.d-3,aexp_frw,hexp_frw,tau_frw,t_frw,n_frw,time_tot)

     ! Find neighboring expansion factors
     i=1
     do while(aexp_frw(i)>aexp.and.i<n_frw)
        i=i+1
     end do
     ! Interploate time
     time_simu=t_frw(i)*(t-tau_frw(i-1))/(tau_frw(i)-tau_frw(i-1))+ &
          & t_frw(i-1)*(t-tau_frw(i))/(tau_frw(i-1)-tau_frw(i))
     write(*,*)'Age simu=',(time_tot+time_simu)/(h0*1d5/3.08d24)/(365.*24.*3600.*1d9)
     time_uni=(time_tot+time_simu)! ADD THIS TO CONVERT TO GYR: /(h0*1d5/3.08d24)/(365.*24.*3600.*1d9)
  else
     time_simu=t
     time_uni=t!*scale_t
  endif

  npart_actual=0
  denspartcount=0
  nstar_actual=0

  if(partok) then
     call readpart(ncpu,ncpu_read,cpu_list,ndim,repository,metal,star,sink,&
          & lmin,lmax,xmin,xmax,ymin,ymax,zmin,zmax,nmin,nmax,npart_actual,&
          & ndm_actual,nstar_actual)
     !NOTE: READPART SELECTS ONLY PARTICLES WITH ID>0.

     write(*,*)'Number of DM particles in the selected box: ', npart_actual-nstar_actual
     write(*,*)'Number of star particles in the selected box: ', nstar_actual

     do i=1,npart_actual
        if((ageout(i).ne.0.0d0))then
           if(cosmo)then
              iii=1 ! Compute star formation time
              do while(tau_frw(iii)>ageout(i).and.iii<n_frw)
                 iii=iii+1
              end do
              time=t_frw(iii)*(ageout(i)-tau_frw(iii-1))/(tau_frw(iii)-tau_frw(iii-1))+ &
                   & t_frw(iii-1)*(ageout(i)-tau_frw(iii))/(tau_frw(iii-1)-tau_frw(iii))
              time=(time_tot+time)!ADD THIS TO CONVERT TO YR: /(h0*1d5/3.08d24)/(365.*24.*3600.)
              ageout(i)=time ! Replace age with formation time
           else
              ageout(i)=ageout(i)!ADD THIS TO CONVERT TO YR: *scale_t/(365.*24.*3600.)
           end if
        end if
     end do

     mdm=msph ! IN THE OLD VERSION: mdm*omega_b/(omega_m-omega_b)
     if(mres.ne.0.d0)mdm=max(mdm,mres)
     if(mdm.eq.0)then
        write(*,*)'Use option -mre, please! STOP!'
        stop
     end if
  else
     mdm=msph ! IN THE OLD VERSION: 1d0/(2d0**(3d0*dble(lmax)))
     if(mres.ne.0.d0)mdm=max(mdm,mres)
     if(mdm.eq.0)then
        write(*,*)'Use option -mre, please! STOP!'
        stop
     end if
  end if

  if(hydrok.and.gas)then
     call gaspart3(ncpu,ncpu_read,cpu_list,repository,ordering,&
          & ndummypart,facdens,levelsel,lmax,xmin,xmax,ymin,&
          & ymax,zmin,zmax,mdm,partmass,averdens,&
          & denspartcount)
     nmin=1
     nmax=denspartcount
  end if

  !-------------------------------------------------------------
  !  Writing output tipsy file in ascii format
  !-------------------------------------------------------------
  if(filetype .EQ. 'ascii')then

  write(*,*)'Outputing data in tipsy ASCII format'

  if(gas)then

     open(66,file=outfich,status='unknown',form='formatted')
     if(do_id)then
        open(55,file='partID_'//TRIM(nchar),status='unknown',form='formatted')
     endif

     three=3

     !HEADER
     write(66,*)npart_actual+denspartcount,denspartcount,nstar_actual
     write(66,*)three
     write(66,*)time_uni

     if(do_id)then
        write(55,*)npart_actual+denspartcount,denspartcount,nstar_actual
        write(55,*)three
        write(55,*)time_uni
     endif

     !MASSES AND IDS
     if(hydrok)then
        do i=1,nmax-nmin+1
           write(66,*)partmass
        end do
     end if
     if(partok)then
        do i=1,npart_actual
           if(.not.star)write(66,*)mout(i)
           if((star.or.sink).and.ageout(i)==0.d0)write(66,*)mout(i)
        enddo
        if(do_id)then
           do i=1,npart_actual
              if(.not.star)write(55,*)idout(i)
              if((star.or.sink).and.ageout(i)==0.d0)write(55,*)idout(i)
           enddo
        endif
        if(star.and.nstar_actual>0)then
           do i=1,npart_actual
              if(ageout(i)/=0.d0)write(66,*)mout(i)
           enddo
           if(do_id)then
              do i=1,npart_actual
                 if(ageout(i)/=0.d0)write(55,*)idout(i)
              enddo
           endif
        endif
     end if

     !X COORDINATE
     if(hydrok)then
        do i=1,nmax-nmin+1
           if(varp(i,1)>=facdens*averdens)write(66,*)xp(i,1)
        end do
     end if
     if(partok)then
        do i=1,npart_actual
           if(.not.star)write(66,*)xout(i,1)
           if((star.or.sink).and.ageout(i)==0.d0)write(66,*)xout(i,1)
        enddo
        if(star.and.nstar_actual>0)then
           do i=1,npart_actual
              if(ageout(i)/=0.d0)write(66,*)xout(i,1)
           enddo
        endif
     endif

     !Y COORDINATE
     if(hydrok)then
        do i=1,nmax-nmin+1
           if(varp(i,1)>=facdens*averdens.and.ndim>=2)write(66,*)xp(i,2)
           if(varp(i,1)>=facdens*averdens.and.ndim<2)write(66,*)(ymin+ymax)/2
        end do
     end if
     if(partok)then
        do i=1,npart_actual
           if((.not.star).and.ndim>=2)write(66,*)xout(i,2)
           if((star.or.sink).and.ageout(i)==0.d0.and.ndim>=2)write(66,*)xout(i,2)
           if((.not.star).and.ndim<2)write(66,*)(ymin+ymax)/2
           if((star.or.sink).and.ageout(i)==0.d0.and.ndim<2)write(66,*)(ymin+ymax)/2
        enddo
        if(star.and.nstar_actual>0)then
           do i=1,npart_actual
              if(ageout(i)/=0.d0.and.ndim>=2)write(66,*)xout(i,2)
              if(ageout(i)/=0.d0.and.ndim<2)write(66,*)(ymin+ymax)/2
           enddo
        endif
     end if

     !Z COORDINATE
     if(hydrok)then
        do i=1,nmax-nmin+1
           if(varp(i,1)>=facdens*averdens.and.ndim>=3)write(66,*)xp(i,3)
           if(varp(i,1)>=facdens*averdens.and.ndim<3)write(66,*)(zmin+zmax)/2
        end do
     end if
     if(partok)then
        do i=1,npart_actual
           if((.not.star).and.ndim>=3)write(66,*)xout(i,3)
           if((star.or.sink).and.ageout(i)==0.d0.and.ndim>=3)write(66,*)xout(i,3)
           if((.not.star).and.ndim<3)write(66,*)(zmin+zmax)/2
           if((star.or.sink).and.ageout(i)==0.d0.and.ndim<3)write(66,*)(zmin+zmax)/2
        enddo
        if(star.and.nstar_actual>0)then
           do i=1,npart_actual
              if(ageout(i)/=0.d0.and.ndim>=3)write(66,*)xout(i,3)
              if(ageout(i)/=0.d0.and.ndim<3)write(66,*)(zmin+zmax)/2
           enddo
        endif
     end if

     !VELOCITY

     dummy=0.d0

     !V_X
     if(hydrok)then
        do i=1,nmax-nmin+1
           if(varp(i,1)>=facdens*averdens)write(66,*)varp(i,2)
        end do
     end if
     if(partok)then
        do i=1,npart_actual
           if(.not.star)write(66,*)vout(i,1)
           if((star.or.sink).and.ageout(i)==0.d0)write(66,*)vout(i,1)
        enddo
        if(star.and.nstar_actual>0)then
           do i=1,npart_actual
              if(ageout(i)/=0.d0)write(66,*)vout(i,1)
           enddo
        endif
     end if

     !V_Y
     if(hydrok)then
        do i=1,nmax-nmin+1
           if(varp(i,1)>=facdens*averdens.and.ndim>=2)write(66,*)varp(i,3)
           if(varp(i,1)>=facdens*averdens.and.ndim<2)write(66,*)dummy
        end do
     end if
     if(partok)then
        do i=1,npart_actual
           if((.not.star).and.ndim>=2)write(66,*)vout(i,2)
           if((star.or.sink).and.ageout(i)==0.d0.and.ndim>=2)write(66,*)vout(i,2)
           if((.not.star).and.ndim<2)write(66,*)dummy
           if((star.or.sink).and.ageout(i)==0.d0.and.ndim<2)write(66,*)dummy
        enddo
        if(star.and.nstar_actual>0)then
           do i=1,npart_actual
              if(ageout(i)/=0.d0.and.ndim>=2)write(66,*)vout(i,2)
              if(ageout(i)/=0.d0.and.ndim<2)write(66,*)dummy
           enddo
        endif
     endif

     !V_Z
     if(hydrok)then
        do i=1,nmax-nmin+1
           if(varp(i,1)>=facdens*averdens.and.ndim>=3)write(66,*)varp(i,4)
           if(varp(i,1)>=facdens*averdens.and.ndim<3)write(66,*)dummy
        end do
     end if
     if(partok)then
        do i=1,npart_actual
           if((.not.star).and.ndim>=3)write(66,*)vout(i,3)
           if((star.or.sink).and.ageout(i)==0.d0.and.ndim>=3)write(66,*)vout(i,3)
           if((.not.star).and.ndim<3)write(66,*)dummy
           if((star.or.sink).and.ageout(i)==0.d0.and.ndim<3)write(66,*)dummy
        enddo
        if(star.and.nstar_actual>0)then
           do i=1,npart_actual
              if(ageout(i)/=0.d0.and.ndim>=3)write(66,*)vout(i,3)
              if(ageout(i)/=0.d0.and.ndim<3)write(66,*)dummy
           enddo
        endif
     endif

     !DUMMY GRAVITATIONAL SOFTENING FOR DARK AND STARS
     if(partok)then
        dummy=boxlen/2**lmax !THIS IS A DUMMY VALUE: IT CORRESPONDS TO THE CELL SIZE AT THE MAXIMUM LEVEL.
        do i=1,ndm_actual
           write(66,*)dummy
        end do
        if(star.and.nstar_actual>0)then
           do i=1,nstar_actual
              write(66,*)dummy
           end do
        end if
     end if

     !GAS DENSITY, TEMPERATURE, DUMMY SPH SMOOTHING LENGTH & GAS METALLICITY.
     if(hydrok)then
        do i=1,nmax-nmin+1
           if(varp(i,1)>=facdens*averdens)write(66,*)varp(i,1)
        end do
        do i=1,nmax-nmin+1
           if(varp(i,1)>=facdens*averdens.and.ndim==3.and.(.not.mhd))write(66,*)varp(i,5)/(gamma-1.d0)/varp(i,1) !THIS IS P/RHO=(k_b*T)/(mu*m_h)
           if(varp(i,1)>=facdens*averdens.and.ndim==2.and.(.not.mhd))write(66,*)varp(i,4)/(gamma-1.d0)/varp(i,1) !THIS IS P/RHO=(k_b*T)/(mu*m_h)
           if(varp(i,1)>=facdens*averdens.and.ndim==1.and.(.not.mhd))write(66,*)varp(i,3)/(gamma-1.d0)/varp(i,1) !THIS IS P/RHO=(k_b*T)/(mu*m_h)
           if(varp(i,1)>=facdens*averdens.and.(mhd))write(66,*)varp(i,11)/(gamma-1.d0)/varp(i,1) !THIS IS P/RHO=(k_b*T)/(mu*m_h)
        end do
        dummy=boxlen/2**lmax !THIS IS A DUMMY VALUE: IT CORRESPONDS TO THE CELL SIZE AT THE MAXIMUM LEVEL.
        do i=1,nmax-nmin+1
           write(66,*)dummy
        end do
        do i=1,nmax-nmin+1
           if(varp(i,1)>=facdens*averdens.and.metgas.and.ndim==3.and.(.not.mhd))write(66,*)varp(i,6)
           if(varp(i,1)>=facdens*averdens.and.metgas.and.ndim==2.and.(.not.mhd))write(66,*)varp(i,5)
           if(varp(i,1)>=facdens*averdens.and.metgas.and.ndim==1.and.(.not.mhd))write(66,*)varp(i,4)
           if(varp(i,1)>=facdens*averdens.and.metgas.and.(mhd))write(66,*)varp(i,12)
        end do
     end if

     if(star.and.nstar_actual>0)then
        if(metal)then
           do i=1,npart_actual
              if(ageout(i)/=0.d0)write(66,*)metout(i)
           end do
        else
           dummy=0.d0
           do i=1,nstar_actual
              write(66,*)dummy
           end do
        end if
        do i=1,npart_actual
           if(ageout(i)/=0.d0)write(66,*)ageout(i)
        end do
     end if

     dummy=1.d0
     do i=1,npart_actual+denspartcount
        write(66,*)dummy
     end do

     close(66)

     if(hydrok)deallocate(xp,varp)

     if(partok)then
        deallocate(xout,vout,mout,idout)
        if(star.and.nstar_actual>0)then
           deallocate(ageout)
           if(metal)deallocate(metout)
        end if
     end if

  else

     open(66,file=outfich,status='unknown',form='formatted')
     if(do_id)then
        open(55,file='partID_'//TRIM(nchar),status='unknown',form='formatted')
     endif

     three=3

     !HEADER
     write(66,*)npart_actual,0,nstar_actual
     write(66,*)three
     write(66,*)time_uni

     if(do_id)then
        write(55,*)npart_actual,0,nstar_actual
        write(55,*)three
        write(55,*)time_uni
     endif

     !MASSES AND IDS
     if(partok)then
        do i=1,npart_actual
           if(.not.star)write(66,*)mout(i)
           if((star.or.sink).and.ageout(i)==0.d0)write(66,*)mout(i)
        enddo
        if(do_id)then
           do i=1,npart_actual
              if(.not.star)write(55,*)idout(i)
              if((star.or.sink).and.ageout(i)==0.d0)write(55,*)idout(i)
           enddo
        endif
        if(star.and.nstar_actual>0)then
           do i=1,npart_actual
              if(ageout(i)/=0.d0)write(66,*)mout(i)
           enddo
           if(do_id)then
              do i=1,npart_actual
                 if(ageout(i)/=0.d0)write(55,*)idout(i)
              enddo
           endif
        endif
     end if

     !X COORDINATE
     if(partok)then
        do i=1,npart_actual
           if(.not.star)write(66,*)xout(i,1)
           if((star.or.sink).and.ageout(i)==0.d0)write(66,*)xout(i,1)
        enddo
        if(star.and.nstar_actual>0)then
           do i=1,npart_actual
              if(ageout(i)/=0.d0)write(66,*)xout(i,1)
           enddo
        endif
     endif

     !Y COORDINATE
     if(partok)then
        do i=1,npart_actual
           if((.not.star).and.ndim>=2)write(66,*)xout(i,2)
           if((star.or.sink).and.ageout(i)==0.d0.and.ndim>=2)write(66,*)xout(i,2)
           if((.not.star).and.ndim<2)write(66,*)(ymin+ymax)/2
           if((star.or.sink).and.ageout(i)==0.d0.and.ndim<2)write(66,*)(ymin+ymax)/2
        enddo
        if(star.and.nstar_actual>0)then
           do i=1,npart_actual
              if(ageout(i)/=0.d0.and.ndim>=2)write(66,*)xout(i,2)
              if(ageout(i)/=0.d0.and.ndim<2)write(66,*)(ymin+ymax)/2
           enddo
        endif
     end if

     !Z COORDINATE
     if(partok)then
        do i=1,npart_actual
           if((.not.star).and.ndim>=3)write(66,*)xout(i,3)
           if((star.or.sink).and.ageout(i)==0.d0.and.ndim>=3)write(66,*)xout(i,3)
           if((.not.star).and.ndim<3)write(66,*)(zmin+zmax)/2
           if((star.or.sink).and.ageout(i)==0.d0.and.ndim<3)write(66,*)(zmin+zmax)/2
        enddo
        if(star.and.nstar_actual>0)then
           do i=1,npart_actual
              if(ageout(i)/=0.d0.and.ndim>=3)write(66,*)xout(i,3)
              if(ageout(i)/=0.d0.and.ndim<3)write(66,*)(zmin+zmax)/2
           enddo
        endif
     end if

     !VELOCITY

     dummy=0.d0

     !V_X
     if(partok)then
        do i=1,npart_actual
           if(.not.star)write(66,*)vout(i,1)
           if((star.or.sink).and.ageout(i)==0.d0)write(66,*)vout(i,1)
        enddo
        if(star.and.nstar_actual>0)then
           do i=1,npart_actual
              if(ageout(i)/=0.d0)write(66,*)vout(i,1)
           enddo
        endif
     end if

     !V_Y
     if(partok)then
        do i=1,npart_actual
           if((.not.star).and.ndim>=2)write(66,*)vout(i,2)
           if((star.or.sink).and.ageout(i)==0.d0.and.ndim>=2)write(66,*)vout(i,2)
           if((.not.star).and.ndim<2)write(66,*)dummy
           if((star.or.sink).and.ageout(i)==0.d0.and.ndim<2)write(66,*)dummy
        enddo
        if(star.and.nstar_actual>0)then
           do i=1,npart_actual
              if(ageout(i)/=0.d0.and.ndim>=2)write(66,*)vout(i,2)
              if(ageout(i)/=0.d0.and.ndim<2)write(66,*)dummy
           enddo
        endif
     endif

     !V_Z
     if(partok)then
        do i=1,npart_actual
           if((.not.star).and.ndim>=3)write(66,*)vout(i,3)
           if((star.or.sink).and.ageout(i)==0.d0.and.ndim>=3)write(66,*)vout(i,3)
           if((.not.star).and.ndim<3)write(66,*)dummy
           if((star.or.sink).and.ageout(i)==0.d0.and.ndim<3)write(66,*)dummy
        enddo
        if(star.and.nstar_actual>0)then
           do i=1,npart_actual
              if(ageout(i)/=0.d0.and.ndim>=3)write(66,*)vout(i,3)
              if(ageout(i)/=0.d0.and.ndim<3)write(66,*)dummy
           enddo
        endif
     endif

     !DUMMY GRAVITATIONAL SOFTENING FOR DARK AND STARS
     if(partok)then
        dummy=boxlen/2**lmax !THIS IS A DUMMY VALUE: IT CORRESPONDS TO THE CELL SIZE AT THE MAXIMUM LEVEL.
        do i=1,ndm_actual
           write(66,*)dummy
        end do
        if(star.and.nstar_actual>0)then
           do i=1,nstar_actual
              write(66,*)dummy
           end do
        end if
     end if

     if(star.and.nstar_actual>0)then
        if(metal)then
           do i=1,npart_actual
              if(ageout(i)/=0.d0)write(66,*)metout(i)
           end do
        else
           dummy=0.d0
           do i=1,nstar_actual
              write(66,*)dummy
           end do
        end if
        do i=1,npart_actual
           if(ageout(i)/=0.d0)write(66,*)ageout(i)
        end do
     end if

     dummy=1.d0
     do i=1,npart_actual
        write(66,*)dummy
     end do

     close(66)

  end if

  !-------------------------------------------------------------
  !  Writing output tipsy file in binary format
  !-------------------------------------------------------------
  else

     write(*,*)'Outputing data in tipsy BINARY format'

     open(66,file=outfich,status='unknown',form='unformatted',access='direct',recl=1)

     three=3

     !HEADER
     write(66,rec=1)real(time_uni)
     write(66,rec=2)real(time_uni)
     write(66,rec=3)npart_actual+denspartcount
     write(66,rec=4)three
     write(66,rec=5)denspartcount
     write(66,rec=6)ndm_actual
     write(66,rec=7)nstar_actual

     write(*,*)'Header done'

     !GAS PARTICLES
     if(hydrok)then
        do i=1,denspartcount
           iskip=8+(i-1)*12
           write(66,rec=iskip+1)real(partmass)
           write(66,rec=iskip+2)real(xp(i,1))
           write(66,rec=iskip+3)real(xp(i,2))
           write(66,rec=iskip+4)real(xp(i,3))
           write(66,rec=iskip+5)real(varp(i,2))
           write(66,rec=iskip+6)real(varp(i,3))
           write(66,rec=iskip+7)real(varp(i,4))
           write(66,rec=iskip+8)real(varp(i,1))
           write(66,rec=iskip+9)real(varp(i,5)/varp(i,1)/(gamma-1.0))
           write(66,rec=iskip+10)real(boxlen/2**lmax)
           if(metal)then
              write(66,rec=iskip+11)real(varp(i,6))
           else
              write(66,rec=iskip+11)real(0.*boxlen)
           endif
           write(66,rec=iskip+12)real(0.*boxlen)
        end do

        write(*,*)'Gas particles done'

     endif

     !DM PARTICLES
     if(partok)then
        if(ndm_actual>0)then
           j=1
           do i=1,npart_actual
              iskip=8+denspartcount*12+(j-1)*9
              if(ageout(i)==0)then
                 write(66,rec=iskip+1)real(mout(i))
                 write(66,rec=iskip+2)real(xout(i,1))
                 write(66,rec=iskip+3)real(xout(i,2))
                 write(66,rec=iskip+4)real(xout(i,3))
                 write(66,rec=iskip+5)real(vout(i,1))
                 write(66,rec=iskip+6)real(vout(i,2))
                 write(66,rec=iskip+7)real(vout(i,3))
                 write(66,rec=iskip+10)real(boxlen/2**lmax)
                 write(66,rec=iskip+11)real(0.*boxlen)
                 j=j+1
              endif
           end do

           write(*,*)'Dark matter particles done'

        endif
     endif

     !STAR PARTICLES
     if(partok)then
        if(nstar_actual>0)then
           j=1
           do i=1,npart_actual
              iskip=8+denspartcount*12+ndm_actual*9+(j-1)*11
              if(ageout(i)/=0)then
                 write(66,rec=iskip+1)real(mout(i))
                 write(66,rec=iskip+2)real(xout(i,1))
                 write(66,rec=iskip+3)real(xout(i,2))
                 write(66,rec=iskip+4)real(xout(i,3))
                 write(66,rec=iskip+5)real(vout(i,1))
                 write(66,rec=iskip+6)real(vout(i,2))
                 write(66,rec=iskip+7)real(vout(i,3))
                 if(metal)then
                    write(66,rec=iskip+11)real(metout(i))
                 else
                    write(66,rec=iskip+11)real(0.*boxlen)
                 endif
                 write(66,rec=iskip+11)real(ageout(i))
                 write(66,rec=iskip+10)real(boxlen/2**lmax)
                 write(66,rec=iskip+11)real(0.*boxlen)
                 j=j+1
              endif
           end do

           write(*,*)'New star particles done'

        endif
     endif

     close(66)

     open(66,file=outfich,status='unknown',form='unformatted',access='direct',recl=2)
     write(66,rec=1)time_uni
     close(66)

  endif

  write(*,*)'File dump completed'

  if(hydrok)deallocate(xp,varp)

  if(partok)then
     deallocate(xout,vout,mout,idout)
     if(star.and.nstar_actual>0)then
        deallocate(ageout)
        if(metal)deallocate(metout)
     end if
  end if

contains

  subroutine read_params

    implicit none

    integer       :: i,n
    integer       :: iargc
    character(len=4)   :: opt
    character(len=128) :: arg
    LOGICAL       :: bad, ok

    n = iargc()
    if (n < 4) then
       print *, 'usage: ramses2tipsy -inp  input_dir'
       print *, '                 -out  output_file'
       print *, '                 [-xmi xmin] '
       print *, '                 [-xma xmax] '
       print *, '                 [-ymi ymin] '
       print *, '                 [-yma ymax] '
       print *, '                 [-zmi zmin] '
       print *, '                 [-zma zmax] '
       print *, '                 [-mre gas particle mass] '
       print *, '                 [-fil filetype=bin or ascii] '
       print *, '                 [-cos cosmo run?] '
       print *, '                 [-pid store particle id?] '
       print *, '                 [-mhd mhd?] '
       print *, '                 [-gas store gas?] '
       print *, 'ex: ramses2tipsy -inp output_00001 -out cube.dat'// &
            &   ' -xmi 0.1 -xma 0.7'
       print *, ' '

       stop
    end if

    do i = 1,n,2
       call getarg(i,opt)
       if (i == n) then
          print '("option ",a2," has no argument")', opt
          stop 2
       end if
       call getarg(i+1,arg)
       select case (opt)
       case ('-inp')
          repository = trim(arg)
       case ('-out')
          outfich = trim(arg)
       case ('-xmi')
          read (arg,*) xmin
       case ('-xma')
          read (arg,*) xmax
       case ('-ymi')
          read (arg,*) ymin
       case ('-yma')
          read (arg,*) ymax
       case ('-zmi')
          read (arg,*) zmin
       case ('-zma')
          read (arg,*) zmax
       case ('-cos')
          read (arg,*) cosmo
       case ('-mhd')
          read (arg,*) mhd
       case ('-fil')
          filetype = trim(arg)
       case ('-pid')
          read (arg,*) do_id
       case ('-gas')
          read (arg,*) gas
       case ('-mre')
          read (arg,*) mres
       case default
          print '("unknown option ",a2," ignored")', opt
       end select
    end do

    return

  end subroutine read_params

end program ramses2tipsy

!================================================================
!================================================================
!================================================================
!================================================================

subroutine friedman(O_mat_0,O_vac_0,O_k_0,alpha,axp_min, &
     & axp_out,hexp_out,tau_out,t_out,ntable,age_tot)

  implicit none
  integer::ntable
  real(kind=8)::O_mat_0, O_vac_0, O_k_0
  real(kind=8)::alpha,axp_min,age_tot
  real(kind=8),dimension(0:ntable)::axp_out,hexp_out,tau_out,t_out
  ! ######################################################!
  ! This subroutine assumes that axp = 1 at z = 0 (today) !
  ! and that t and tau = 0 at z = 0 (today).              !
  ! axp is the expansion factor, hexp the Hubble constant !
  ! defined as hexp=1/axp*daxp/dtau, tau the conformal    !
  ! time, and t the look-back time, both in unit of 1/H0. !
  ! alpha is the required accuracy and axp_min is the     !
  ! starting expansion factor of the look-up table.       !
  ! ntable is the required size of the look-up table.     !
  ! ######################################################!
  real(kind=8)::axp_tau, axp_t
  real(kind=8)::axp_tau_pre, axp_t_pre
  real(kind=8)::dadtau, dadt
  real(kind=8)::dtau,dt
  real(kind=8)::tau,t
  integer::nstep,nout,nskip

!  if( (O_mat_0+O_vac_0+O_k_0) .ne. 1.0D0 )then
!     write(*,*)'Error: non-physical cosmological constants'
!     write(*,*)'O_mat_0,O_vac_0,O_k_0=',O_mat_0,O_vac_0,O_k_0
!     write(*,*)'The sum must be equal to 1.0, but '
!     write(*,*)'O_mat_0+O_vac_0+O_k_0=',O_mat_0+O_vac_0+O_k_0
!     stop
!  end if

  axp_tau = 1.0D0
  axp_t = 1.0D0
  tau = 0.0D0
  t = 0.0D0
  nstep = 0

  do while ( (axp_tau .ge. axp_min) .or. (axp_t .ge. axp_min) )

     nstep = nstep + 1
     dtau = alpha * axp_tau / dadtau(axp_tau,O_mat_0,O_vac_0,O_k_0)
     axp_tau_pre = axp_tau - dadtau(axp_tau,O_mat_0,O_vac_0,O_k_0)*dtau/2.d0
     axp_tau = axp_tau - dadtau(axp_tau_pre,O_mat_0,O_vac_0,O_k_0)*dtau
     tau = tau - dtau

     dt = alpha * axp_t / dadt(axp_t,O_mat_0,O_vac_0,O_k_0)
     axp_t_pre = axp_t - dadt(axp_t,O_mat_0,O_vac_0,O_k_0)*dt/2.d0
     axp_t = axp_t - dadt(axp_t_pre,O_mat_0,O_vac_0,O_k_0)*dt
     t = t - dt

  end do

  age_tot=-t
  write(*,666)-t
  666 format(' Age of the Universe (in unit of 1/H0)=',1pe10.3)

  nskip=nstep/ntable

  axp_t = 1.d0
  t = 0.d0
  axp_tau = 1.d0
  tau = 0.d0
  nstep = 0
  nout=0
  t_out(nout)=t
  tau_out(nout)=tau
  axp_out(nout)=axp_tau
  hexp_out(nout)=dadtau(axp_tau,O_mat_0,O_vac_0,O_k_0)/axp_tau

  do while ( (axp_tau .ge. axp_min) .or. (axp_t .ge. axp_min) )

     nstep = nstep + 1
     dtau = alpha * axp_tau / dadtau(axp_tau,O_mat_0,O_vac_0,O_k_0)
     axp_tau_pre = axp_tau - dadtau(axp_tau,O_mat_0,O_vac_0,O_k_0)*dtau/2.d0
     axp_tau = axp_tau - dadtau(axp_tau_pre,O_mat_0,O_vac_0,O_k_0)*dtau
     tau = tau - dtau

     dt = alpha * axp_t / dadt(axp_t,O_mat_0,O_vac_0,O_k_0)
     axp_t_pre = axp_t - dadt(axp_t,O_mat_0,O_vac_0,O_k_0)*dt/2.d0
     axp_t = axp_t - dadt(axp_t_pre,O_mat_0,O_vac_0,O_k_0)*dt
     t = t - dt

     if(mod(nstep,nskip)==0)then
        nout=nout+1
        t_out(nout)=t
        tau_out(nout)=tau
        axp_out(nout)=axp_tau
        hexp_out(nout)=dadtau(axp_tau,O_mat_0,O_vac_0,O_k_0)/axp_tau
     end if
  end do
  t_out(ntable)=t
  tau_out(ntable)=tau
  axp_out(ntable)=axp_tau
  hexp_out(ntable)=dadtau(axp_tau,O_mat_0,O_vac_0,O_k_0)/axp_tau

end subroutine friedman

function dadtau(axp_tau,O_mat_0,O_vac_0,O_k_0)
  real(kind=8)::dadtau,axp_tau,O_mat_0,O_vac_0,O_k_0
  dadtau = axp_tau*axp_tau*axp_tau *  &
       &   ( O_mat_0 + &
       &     O_vac_0 * axp_tau*axp_tau*axp_tau + &
       &     O_k_0   * axp_tau )
  dadtau = sqrt(dadtau)
  return
end function dadtau

function dadt(axp_t,O_mat_0,O_vac_0,O_k_0)
  real(kind=8)::dadt,axp_t,O_mat_0,O_vac_0,O_k_0
  dadt   = (1.0D0/axp_t)* &
       &   ( O_mat_0 + &
       &     O_vac_0 * axp_t*axp_t*axp_t + &
       &     O_k_0   * axp_t )
  dadt = sqrt(dadt)
  return
end function dadt

max_forks_repo_name: rknop/amuse