1. Matthew Turk
  2. enzo-dengo

Commits

Daniel Reynolds  committed f78b5a3

Updated fortran precision handling inside inits to match new Enzo infrastructure (I forgot to do that before the 2.2 release).

  • Participants
  • Parent commits d992797
  • Branches week-of-code

Comments (0)

Files changed (82)

File src/inits/Makefile

View file
  • Ignore whitespace
 # Implicit rules
 #-----------------------------------------------------------------------
 
-.SUFFIXES: .c .C .src .src90 .f .f90
+.SUFFIXES: .c .C .src .src90 .f .f90 .F .F90
+
+.F.o:
+	@echo "Compiling $<"
+	@rm -f $@
+	@(if [ $(VERBOSE) -eq 0 ]; then \
+	  $(FC) -c -o $@ $(FFLAGS) $(DEFINES) $*.F >& $(OUTPUT) ; \
+	  if [ ! -e $@ ]; then \
+             echo; \
+             echo "$(FC) -c -o $@ $(FFLAGS) $(DEFINES) $*.F"; \
+             echo; \
+             $(FC) -c -o $@ $(FFLAGS) $(DEFINES) $*.F; \
+             echo; \
+             exit 1; \
+          fi ; \
+	else \
+	  $(FC) -c -o $@ $(FFLAGS)  $(DEFINES) $*.F >> $(OUTPUT) 2>&1 ; \
+	  if [ ! -e $@ ]; then \
+	     echo "See $(OUTPUT) for error messages"; \
+	     exit 1; \
+	  fi ; \
+	fi)
+
+.F90.o:
+	@rm -f $@
+	@echo "Compiling $<"
+	@(if [ $(VERBOSE) -eq 0 ]; then \
+	  $(F90) -c -o $@ $(F90FLAGS) $(DEFINES) $*.F90 >& $(OUTPUT) ; \
+	  if [ ! -e $@ ]; then \
+             echo; \
+             echo "$(F90) -c -o $@ $(DEFINES) $(F90FLAGS) $*.F90"; \
+             echo; \
+             $(F90) -c -o $@ $(F90FLAGS) $*.F90; \
+             echo; \
+             exit 1; \
+	  fi ; \
+	else \
+	  $(F90) -c -o $@ $(F90FLAGS) $(DEFINES) $*.F90 >> $(OUTPUT) 2>&1 ; \
+	  if [ ! -e $@ ]; then \
+	     echo "See $(OUTPUT) for error messages"; \
+	     exit 1; \
+	  fi ; \
+	fi)
 
 .src.f:
 	@echo "Compiling $<"

File src/inits/Rmake_field.F

View file
  • Ignore whitespace
+#include "../enzo/fortran.def"
+!=======================================================================
+!//////////////////////  SUBROUTINE MAKE_FIELD  \\\\\\\\\\\\\\\\\\\\\\\\
+
+      subroutine make_field(field, nx, ny, nz, nxmax, nymax, nzmax,
+     &                      in, jn, kn, itype, iseed, box,
+     &                      PSTable, PSMin, PSStep, kfcutoff)
+
+!  COMPUTES RANDOM GAUSSIAN FIELD FROM SPECIFIED POWER SPECTRUM
+!
+!  written by: Greg Bryan
+!  date:       June, 1997
+!  modified:   Robert Harkness
+!  date:       November, 2003
+!
+!  PURPOSE: 
+!
+!  INPUTS:
+!        i,j,kn      = real dimensions of green
+!        nx,ny,nz    = active dimensions of green
+!        nx,y,zmax   = dimensions of k field (for random number countinf)
+!        itype       = field type (0 - density, 1/2/3 - x/y/z displacement)
+!        iseed       = random number seed (negative)
+!        box         = size
+!        PSTable     = Table of precomputed PS values
+!        PSMin       = minimum x value in PSTable
+!        PSStep      = x step in PSTable
+!        kfcutoff    = high k filter (sharp) in units of the fundamental
+!
+!  Outputs:
+!        field       = gaussian random field
+!
+!  LOCALS:
+!        num_dim     = number of dimensions to be used for force law
+!        nx,y,zmid   = midpoint (+1) of the grid in each axis
+!        nx,y,zd2    = number of grid points divided by 2 for each axis
+
+      implicit NONE
+#include "../enzo/fortran_types.def"
+
+!     Arguments
+
+      INTG_PREC :: in, jn, kn, nx, ny, nz, nxmax, nymax, nzmax, 
+     &           itype, iseed, kfcutoff
+      R_PREC ::    field(in, jn, kn), box, 
+     &           PSMin, PSPart, PSStep, PSTable(*)
+
+!     Locals
+
+      INTG_PREC :: i, ii, j, jj, k, kk, index, nxmid, nymid, nzmid
+      R_PREC ::    ang, amp, d3k, dummy, kmodsq, kx, ky, kz, kdir,
+     &           klog, psval, twopi, kcutoffsq, dkx, dky, dkz
+      CMPLX_PREC :: z
+
+      R_PREC ::    ranf_min
+      INTG_PREC :: long_seed
+
+!     External function
+
+      R_PREC ::    ran1, enzo_ranf
+
+!  Define table lookup function
+
+      R_PREC ::    Table1, Table2, Step, Min, Tablex, TableLookUp
+      INTG_PREC :: Tablei
+
+      TableLookUp(Table1, Table2, Step, Min, Tablei, Tablex) = 
+     &            Table1 + (Tablex - REAL(Tablei-1,RKIND)*Step - Min) 
+     &            / Step * (Table2 - Table1)
+
+
+!     Set constants
+
+      twopi  = 8.0_RKIND*atan(1.0_RKIND)
+      dkx    = twopi/box
+      dky    = twopi/(ny*box/nx)
+      dkz    = twopi/(nz*box/nx)
+      d3k    = (twopi/box)**3
+      kcutoffsq = 1.0e30_RKIND
+
+      if (kfcutoff .gt. 0) kcutoffsq = (kfcutoff*dkx)**2
+
+!     Initialize random # generator with random seed
+
+      long_seed = iseed
+
+!     dummy = ran1(iseed)
+
+      call enzo_seed(long_seed)
+
+!     Set maximum wavenumbers in the positive direction
+
+      nxmid  = max(nx/2_IKIND + 1_IKIND, 1_IKIND)
+      nymid  = max(ny/2_IKIND + 1_IKIND, 1_IKIND)
+      nzmid  = max(nz/2_IKIND + 1_IKIND, 1_IKIND)
+
+!     Loop over k
+
+      do k = 1, nzmax
+         kk = k-1
+         if (k .gt. nzmax/2+1) kk = kk - nzmax
+         kz    = REAL(kk,RKIND)*dkz
+
+!        Loop over j
+
+         do j = 1, nymax
+            jj = j-1
+            if (j .gt. nymax/2+1) jj = jj - nymax
+            ky    = REAL(jj,RKIND)*dky
+ 
+!           If this j corresponds to a wavenumber vector that is not
+!           inside field then just skip over it (& eat the random nums).
+
+            if (jj .ge. nymid .or. jj .lt. -ny/2+1 .or.
+     &          kk .ge. nzmid .or. kk .lt. -nz/2+1     ) then
+
+               do i=1,nxmax/2+1
+                  dummy = enzo_ranf()
+                  dummy = enzo_ranf()
+               enddo
+
+               goto 100
+
+            endif
+
+!           Loop over i
+
+            do i = 1, nxmid
+               ii    = i-1
+               kx    = REAL(ii,RKIND)*dkx
+
+!              Compute kmod and lookup value in table
+!              (and convert from power per mode).
+
+               kmodsq  = (kx**2 + ky**2 + kz**2)
+               if (kmodsq .eq. 0) kmodsq = 1.0_RKIND
+               klog   = 0.5_RKIND*log(kmodsq)
+               index = int((klog - PSMin)/PSStep,IKIND)
+               psval = TableLookUp(PSTable(index), PSTable(index+1),
+     &                             PSStep, PSMin, index, klog)
+               psval = psval * d3k
+               if (kmodsq .gt. kcutoffsq) psval = 0.0_RKIND
+
+!              Generate a complex number with random phase and amplitude
+!              Gaussian distributed with a mean of sqrt(psval) with the
+!              Box-Muller method.  Note we have supressed a factor of
+!              sqrt(2) since we must also divide by this factor to account
+!              for the dreary fact that we are really generating two random
+!              fields (if we were doing a complex-to-complex transform
+!              this would show up when we discarded the perfectly
+!              good imaginary component of the transformed field).  whew.
+
+               ranf_min = 1.e-37_RKIND
+
+               ang = twopi*enzo_ranf()
+               amp = sqrt(-log(max(enzo_ranf(),ranf_min)) * psval)
+               z   = cmplx(cos(ang), sin(ang)) * amp
+
+!              Process this on the basis of itype:
+!                0)   density field - just leave it be.
+!                1-3) displacement field - multiply by vec(k)/k^2
+!                     (and then convert from Mpc to fraction of box).
+
+               if (itype .ne. 0) then
+                  if (itype .eq. 1) kdir = kx
+                  if (itype .eq. 2) kdir = ky
+                  if (itype .eq. 3) kdir = kz
+                  z = z * cmplx(0._RKIND,1._RKIND) * kdir
+     &                 / (kmodsq**2 * box)
+               endif
+
+!              Set the complex field
+
+               field(i*2-1,j,k) = REAL(z,RKIND)
+               field(i*2  ,j,k) = IMAG(z)
+
+            enddo
+
+!           Now loop over the rest of the kx space to use up the
+!           allotted number of random numbers
+
+            do i=nxmid+1, nxmax/2+1
+               dummy = enzo_ranf()
+               dummy = enzo_ranf()
+            enddo
+
+ 100        continue
+
+         enddo
+      enddo
+
+!     Clear the zero wavenumber position
+
+      field(1,1,1) = 0.0_RKIND
+      field(2,1,1) = 0.0_RKIND
+
+!     Adjust the field to satisfy that conjugate relations that
+!     are implied by a zero imaginary part.
+
+      call adjfft(field, nx, ny, nz, in, jn)
+
+      return
+      end

File src/inits/Rmake_field.src

  • Ignore whitespace
-!=======================================================================
-!//////////////////////  SUBROUTINE MAKE_FIELD  \\\\\\\\\\\\\\\\\\\\\\\\
-
-      subroutine make_field(field, nx, ny, nz, nxmax, nymax, nzmax,
-     &                      in, jn, kn, itype, iseed, box,
-     &                      PSTable, PSMin, PSStep, kfcutoff)
-
-!  COMPUTES RANDOM GAUSSIAN FIELD FROM SPECIFIED POWER SPECTRUM
-!
-!  written by: Greg Bryan
-!  date:       June, 1997
-!  modified:   Robert Harkness
-!  date:       November, 2003
-!
-!  PURPOSE: 
-!
-!  INPUTS:
-!        i,j,kn      = real dimensions of green
-!        nx,ny,nz    = active dimensions of green
-!        nx,y,zmax   = dimensions of k field (for random number countinf)
-!        itype       = field type (0 - density, 1/2/3 - x/y/z displacement)
-!        iseed       = random number seed (negative)
-!        box         = size
-!        PSTable     = Table of precomputed PS values
-!        PSMin       = minimum x value in PSTable
-!        PSStep      = x step in PSTable
-!        kfcutoff    = high k filter (sharp) in units of the fundamental
-!
-!  Outputs:
-!        field       = gaussian random field
-!
-!  LOCALS:
-!        num_dim     = number of dimensions to be used for force law
-!        nx,y,zmid   = midpoint (+1) of the grid in each axis
-!        nx,y,zd2    = number of grid points divided by 2 for each axis
-
-      implicit NONE
-
-!     Arguments
-
-      integer :: in, jn, kn, nx, ny, nz, nxmax, nymax, nzmax, 
-     &           itype, iseed, kfcutoff
-      real ::    field(in, jn, kn), box, 
-     &           PSMin, PSPart, PSStep, PSTable(*)
-
-!     Locals
-
-      integer :: i, ii, j, jj, k, kk, index, nxmid, nymid, nzmid
-      real ::    ang, amp, d3k, dummy, kmodsq, kx, ky, kz, kdir,
-     &           klog, psval, twopi, kcutoffsq, dkx, dky, dkz
-      complex :: z
-
-      real ::    ranf_min
-      integer :: long_seed
-
-!     External function
-
-      real ::    ran1, enzo_ranf
-
-!  Define table lookup function
-
-      real ::    Table1, Table2, Step, Min, Tablex, TableLookUp
-      integer :: Tablei
-
-      TableLookUp(Table1, Table2, Step, Min, Tablei, Tablex) = 
-     &            Table1 + (Tablex - real(Tablei-1)*Step - Min) 
-     &            / Step * (Table2 - Table1)
-
-
-!     Set constants
-
-      twopi  = 8.0*atan(1.0)
-      dkx    = twopi/box
-      dky    = twopi/(ny*box/nx)
-      dkz    = twopi/(nz*box/nx)
-      d3k    = (twopi/box)**3
-      kcutoffsq = 1.0e+30
-
-      if (kfcutoff .gt. 0) kcutoffsq = (kfcutoff*dkx)**2
-
-!     Initialize random # generator with random seed
-
-      long_seed = iseed
-
-!     dummy = ran1(iseed)
-
-      call enzo_seed(long_seed)
-
-!     Set maximum wavenumbers in the positive direction
-
-      nxmid  = max(nx/2 + 1, 1)
-      nymid  = max(ny/2 + 1, 1)
-      nzmid  = max(nz/2 + 1, 1)
-
-!     Loop over k
-
-      do k = 1, nzmax
-         kk = k-1
-         if (k .gt. nzmax/2+1) kk = kk - nzmax
-         kz    = real(kk)*dkz
-
-!        Loop over j
-
-         do j = 1, nymax
-            jj = j-1
-            if (j .gt. nymax/2+1) jj = jj - nymax
-            ky    = real(jj)*dky
- 
-!           If this j corresponds to a wavenumber vector that is not
-!           inside field then just skip over it (& eat the random nums).
-
-            if (jj .ge. nymid .or. jj .lt. -ny/2+1 .or.
-     &          kk .ge. nzmid .or. kk .lt. -nz/2+1     ) then
-
-               do i=1,nxmax/2+1
-                  dummy = enzo_ranf()
-                  dummy = enzo_ranf()
-               enddo
-
-               goto 100
-
-            endif
-
-!           Loop over i
-
-            do i = 1, nxmid
-               ii    = i-1
-               kx    = real(ii)*dkx
-
-!              Compute kmod and lookup value in table
-!              (and convert from power per mode).
-
-               kmodsq  = (kx**2 + ky**2 + kz**2)
-               if (kmodsq .eq. 0) kmodsq = 1.0
-               klog   = 0.5*log(kmodsq)
-               index = int((klog - PSMin)/PSStep)
-               psval = TableLookUp(PSTable(index), PSTable(index+1),
-     &                             PSStep, PSMin, index, klog)
-               psval = psval * d3k
-               if (kmodsq .gt. kcutoffsq) psval = 0.0
-
-!              Generate a complex number with random phase and amplitude
-!              Gaussian distributed with a mean of sqrt(psval) with the
-!              Box-Muller method.  Note we have supressed a factor of
-!              sqrt(2) since we must also divide by this factor to account
-!              for the dreary fact that we are really generating two random
-!              fields (if we were doing a complex-to-complex transform
-!              this would show up when we discarded the perfectly
-!              good imaginary component of the transformed field).  whew.
-
-               ranf_min = 1.0e-37
-
-               ang = twopi*enzo_ranf()
-               amp = sqrt(-log(max(enzo_ranf(),ranf_min)) * psval)
-               z   = cmplx(cos(ang), sin(ang)) * amp
-
-!              Process this on the basis of itype:
-!                0)   density field - just leave it be.
-!                1-3) displacement field - multiply by vec(k)/k^2
-!                     (and then convert from Mpc to fraction of box).
-
-               if (itype .ne. 0) then
-                  if (itype .eq. 1) kdir = kx
-                  if (itype .eq. 2) kdir = ky
-                  if (itype .eq. 3) kdir = kz
-                  z = z * cmplx(0.0,1.0) * kdir / (kmodsq**2 * box)
-               endif
-
-!              Set the complex field
-
-               field(i*2-1,j,k) = real(z)
-               field(i*2  ,j,k) = imag(z)
-
-            enddo
-
-!           Now loop over the rest of the kx space to use up the
-!           allotted number of random numbers
-
-            do i=nxmid+1, nxmax/2+1
-               dummy = enzo_ranf()
-               dummy = enzo_ranf()
-            enddo
-
- 100        continue
-
-         enddo
-      enddo
-
-!     Clear the zero wavenumber position
-
-      field(1,1,1) = 0.0
-      field(2,1,1) = 0.0
-
-!     Adjust the field to satisfy that conjugate relations that
-!     are implied by a zero imaginary part.
-
-      call adjfft(field, nx, ny, nz, in, jn)
-
-      return
-      end

File src/inits/Rmake_field_kpreserving.F

View file
  • Ignore whitespace
+#include "../enzo/fortran.def"
+!=======================================================================
+!//////////////////////  SUBROUTINE MAKE_FIELD  \\\\\\\\\\\\\\\\\\\\\\\\
+
+      subroutine make_field_kpreserving(field, nx, ny, nz, 
+     &                      in, jn, kn, itype, iseed, box,
+     &                      PSTable, PSMin, PSStep, kfcutoff)
+
+!  COMPUTES RANDOM GAUSSIAN FIELD FROM SPECIFIED POWER SPECTRUM
+!
+!  written by: Greg Bryan
+!  date:       June, 1997
+!  modified:   Robert Harkness
+!  date:       November, 2003
+!
+!  PURPOSE: 
+!
+!  INPUTS:
+!        i,j,kn      = real dimensions of green
+!        nx,ny,nz    = active dimensions of green
+!        itype       = field type (0 - density, 1/2/3 - x/y/z displacement)
+!        iseed       = random number seed (negative)
+!        box         = size
+!        PSTable     = Table of precomputed PS values
+!        PSMin       = minimum x value in PSTable
+!        PSStep      = x step in PSTable
+!        kfcutoff    = high k filter (sharp) in units of the fundamental
+!
+!  Outputs:
+!        field       = gaussian random field
+!
+!  LOCALS:
+!        num_dim     = number of dimensions to be used for force law
+!        nx,y,zmid   = midpoint (+1) of the grid in each axis
+!        nx,y,zd2    = number of grid points divided by 2 for each axis
+
+      implicit NONE
+#include "../enzo/fortran_types.def"
+
+!     Arguments
+
+      INTG_PREC :: in, jn, kn, nx, ny, nz, nxmax, nymax, nzmax, 
+     &           itype, iseed, kfcutoff
+      R_PREC ::    field(in, jn, kn), box, 
+     &           PSMin, PSPart, PSStep, PSTable(1)
+
+!     Locals
+
+      INTG_PREC :: i, i1, j, j1, n, n1
+      R_PREC ::    dummy, twopi, kcutoffsq, dk
+      CMPLX_PREC :: z
+
+      INTG_PREC :: long_seed
+
+!     External functions
+
+      R_PREC ::    ran1
+
+!     Set constants
+
+      twopi  = 8.0_RKIND*atan(1.0_RKIND)
+      dk     = twopi/box
+      kcutoffsq = 1.0e30_RKIND
+
+      if (kfcutoff .gt. 0) kcutoffsq = (kfcutoff*dk)**2
+
+!     Initialize random # generator with random seed
+
+      long_seed = iseed
+      n = iseed
+!     dummy = ran1(n)
+      call enzo_seed(long_seed)
+
+!     Loop over k-box sizes, so that we fill k-space from low-k outwards
+
+      do n=1,nx/2
+
+         do i=-n+1, n
+            do j=-n+1, n
+
+               i1 = mod(i+nx,nx)+1
+               j1 = mod(j+nx,nx)+1
+               n1 = mod(1_IKIND-n+nx,nx)+1
+
+!              1) +i plane
+
+               call processk(n,i,j, dk, PSMin, PSStep, PSTable, 
+     &                       itype, z, kcutoffsq, box)
+
+               field((n+1)*2-1,i1,j1) = REAL(z,RKIND)
+               field((n+1)*2  ,i1,j1) = imag(z)
+
+!              2) +j and -j plane
+!                 (the i .ne. n is to avoid overlapping with (1))
+
+               if (i .ge. 0 .and. i .ne. n) then
+
+                  call processk(i,n,j, dk, PSMin, PSStep, PSTable, 
+     &                          itype, z, kcutoffsq, box)
+
+                  field(i1*2-1,n+1,j1) = REAL(z,RKIND)
+                  field(i1*2  ,n+1,j1) = imag(z)
+
+                  call processk(i,1_IKIND-n,j, dk, PSMin, PSStep, 
+     &                          PSTable, itype, z, kcutoffsq, box)
+
+                  field(i1*2-1,n1,j1) = REAL(z,RKIND)
+                  field(i1*2  ,n1,j1) = imag(z)
+
+               endif
+
+!              3) +k and -k plane
+!                 (the logic involving j is to avoid overlapping with (2))
+
+               if (i .ge. 0 .and. i .ne. n .and. 
+     &             j .ne. -n+1 .and. j .ne. n) then
+
+                  call processk(i,j,n, dk, PSMin, PSStep, PSTable,
+     &                          itype, z, kcutoffsq, box)
+
+                  field(i1*2-1,j1,n+1) = REAL(z,RKIND)
+                  field(i1*2  ,j1,n+1) = imag(z)
+
+                  call processk(i,j,1_IKIND-n, dk, PSMin, PSStep, 
+     &                          PSTable, itype, z, kcutoffsq, box)
+
+                  field(i1*2-1,j1,n1) = REAL(z,RKIND)
+                  field(i1*2  ,j1,n1) = imag(z)
+
+               endif
+
+            enddo
+         enddo
+
+      enddo
+
+      do i=1, in
+         do j=1, jn
+            do n=1, kn
+               field(i,j,n) = field(i,j,n) * REAL(nx*ny*nz,RKIND)
+            enddo
+         enddo
+      enddo
+
+!     Clear the zero wavenumber position
+
+      field(1,1,1) = 0.0_RKIND
+      field(2,1,1) = 0.0_RKIND
+
+!     Adjust the field to satisfy the conjugate relations that
+!     are implied by a zero imaginary part.
+
+      call adjfft(field, nx, ny, nz, in, jn)
+
+      return
+      end
+
+
+c===================================================================
+
+      subroutine processk(i, j, k, dk, PSMin, PSStep, PSTable, 
+     &                    itype, z, kcutoffsq, box)
+
+      implicit none
+#include "../enzo/fortran_types.def"
+
+!     Parameter
+
+      R_PREC, parameter :: twopi = 2._RKIND*3.14159265358979324_RKIND
+
+!     Arguments
+
+      INTG_PREC :: i, j, k, itype
+      R_PREC ::    dk, PSMin, PSStep, PSTable(*)
+      R_PREC ::    kcutoffsq, box
+      CMPLX_PREC :: z
+
+!     Locals
+
+      INTG_PREC :: index
+      R_PREC :: psval, kdir, klog, ang, amp, kmodsq
+      R_PREC :: ranf_min
+
+!     External function
+
+      R_PREC :: ran1
+      R_PREC :: enzo_ranf
+
+
+
+!     Define table lookup function
+ 
+      R_PREC ::    Table1, Table2, Step, Min, Tablex, TableLookUp
+      INTG_PREC :: Tablei
+
+      TableLookUp(Table1, Table2, Step, Min, Tablei, Tablex) = 
+     &            Table1 + (Tablex - REAL(Tablei-1,RKIND)*Step - Min) 
+     &            / Step * (Table2 - Table1)
+
+
+      kmodsq = max(i**2 + j**2 + k**2, 1)*dk**2
+      klog   = 0.5_RKIND*log(kmodsq)
+      index = int((klog - PSMin)/PSStep)
+      psval = TableLookUp(PSTable(index), PSTable(index+1),
+     &                    PSStep, PSMin, index, klog)
+      psval = psval * dk**3
+
+      if (kmodsq .gt. kcutoffsq) psval = 0.0_RKIND
+
+!     Generate a complex number with random phase and amplitude
+!     Gaussian distributed with a mean of sqrt(psval) with the
+!     Box-Muller method.  Note we have supressed a factor of
+!     sqrt(2) since we must also divide by this factor to account
+!     for the dreary fact that we are really generating two random
+!     fields (if we were doing a complex-to-complex transform
+!     this would show up when we discarded the perfectly
+!     good imaginary component of the transformed field).  whew.
+
+      ranf_min = 1.e-37_RKIND
+
+      ang = twopi*enzo_ranf()
+      amp = sqrt(-log(max(enzo_ranf(),ranf_min)) * psval)
+      z   = cmplx(cos(ang), sin(ang)) * amp
+
+!     Process this on the basis of itype:
+!      0)   density field - just leave it be.
+!      1-3) displacement field - multiply by vec(k)/k^2
+!           (and then convert from Mpc to fraction of box).
+
+      if (itype .ne. 0) then
+         if (itype .eq. 1) kdir = REAL(i,RKIND)*dk
+         if (itype .eq. 2) kdir = REAL(j,RKIND)*dk
+         if (itype .eq. 3) kdir = REAL(k,RKIND)*dk
+         z = z * cmplx(0._RKIND,1._RKIND) * kdir / (kmodsq * box)
+      endif
+
+      return
+      end

File src/inits/Rmake_field_kpreserving.src

  • Ignore whitespace
-!=======================================================================
-!//////////////////////  SUBROUTINE MAKE_FIELD  \\\\\\\\\\\\\\\\\\\\\\\\
-
-      subroutine make_field_kpreserving(field, nx, ny, nz, 
-     &                      in, jn, kn, itype, iseed, box,
-     &                      PSTable, PSMin, PSStep, kfcutoff)
-
-!  COMPUTES RANDOM GAUSSIAN FIELD FROM SPECIFIED POWER SPECTRUM
-!
-!  written by: Greg Bryan
-!  date:       June, 1997
-!  modified:   Robert Harkness
-!  date:       November, 2003
-!
-!  PURPOSE: 
-!
-!  INPUTS:
-!        i,j,kn      = real dimensions of green
-!        nx,ny,nz    = active dimensions of green
-!        itype       = field type (0 - density, 1/2/3 - x/y/z displacement)
-!        iseed       = random number seed (negative)
-!        box         = size
-!        PSTable     = Table of precomputed PS values
-!        PSMin       = minimum x value in PSTable
-!        PSStep      = x step in PSTable
-!        kfcutoff    = high k filter (sharp) in units of the fundamental
-!
-!  Outputs:
-!        field       = gaussian random field
-!
-!  LOCALS:
-!        num_dim     = number of dimensions to be used for force law
-!        nx,y,zmid   = midpoint (+1) of the grid in each axis
-!        nx,y,zd2    = number of grid points divided by 2 for each axis
-
-      implicit NONE
-
-!     Arguments
-
-      integer :: in, jn, kn, nx, ny, nz, nxmax, nymax, nzmax, 
-     &           itype, iseed, kfcutoff
-      real ::    field(in, jn, kn), box, 
-     &           PSMin, PSPart, PSStep, PSTable(1)
-
-!     Locals
-
-      integer :: i, i1, j, j1, n, n1
-      real ::    dummy, twopi, kcutoffsq, dk
-      complex :: z
-
-      integer :: long_seed
-
-!     External functions
-
-      real ::    ran1
-
-!     Set constants
-
-      twopi  = 8.0*atan(1.0)
-      dk     = twopi/box
-      kcutoffsq = 1.0e30
-
-      if (kfcutoff .gt. 0) kcutoffsq = (kfcutoff*dk)**2
-
-!     Initialize random # generator with random seed
-
-      long_seed = iseed
-      n = iseed
-!     dummy = ran1(n)
-      call enzo_seed(long_seed)
-
-!     Loop over k-box sizes, so that we fill k-space from low-k outwards
-
-      do n=1,nx/2
-
-         do i=-n+1, n
-            do j=-n+1, n
-
-               i1 = mod(i+nx,nx)+1
-               j1 = mod(j+nx,nx)+1
-               n1 = mod(1-n+nx,nx)+1
-
-!              1) +i plane
-
-               call processk(n,i,j, dk, PSMin, PSStep, PSTable, 
-     &                       itype, z, kcutoffsq, box)
-
-               field((n+1)*2-1,i1,j1) = real(z)
-               field((n+1)*2  ,i1,j1) = imag(z)
-
-!              2) +j and -j plane
-!                 (the i .ne. n is to avoid overlapping with (1))
-
-               if (i .ge. 0 .and. i .ne. n) then
-
-                  call processk(i,n,j, dk, PSMin, PSStep, PSTable, 
-     &                          itype, z, kcutoffsq, box)
-
-                  field(i1*2-1,n+1,j1) = real(z)
-                  field(i1*2  ,n+1,j1) = imag(z)
-
-                  call processk(i,1-n,j, dk, PSMin, PSStep, PSTable,
-     &                          itype, z, kcutoffsq, box)
-
-                  field(i1*2-1,n1,j1) = real(z)
-                  field(i1*2  ,n1,j1) = imag(z)
-
-               endif
-
-!              3) +k and -k plane
-!                 (the logic involving j is to avoid overlapping with (2))
-
-               if (i .ge. 0 .and. i .ne. n .and. 
-     &             j .ne. -n+1 .and. j .ne. n) then
-
-                  call processk(i,j,n, dk, PSMin, PSStep, PSTable,
-     &                          itype, z, kcutoffsq, box)
-
-                  field(i1*2-1,j1,n+1) = real(z)
-                  field(i1*2  ,j1,n+1) = imag(z)
-
-                  call processk(i,j,1-n, dk, PSMin, PSStep, PSTable,
-     &                          itype, z, kcutoffsq, box)
-
-                  field(i1*2-1,j1,n1) = real(z)
-                  field(i1*2  ,j1,n1) = imag(z)
-
-               endif
-
-            enddo
-         enddo
-
-      enddo
-
-      do i=1, in
-         do j=1, jn
-            do n=1, kn
-               field(i,j,n) = field(i,j,n) * real(nx*ny*nz)
-            enddo
-         enddo
-      enddo
-
-!     Clear the zero wavenumber position
-
-      field(1,1,1) = 0.0
-      field(2,1,1) = 0.0
-
-!     Adjust the field to satisfy the conjugate relations that
-!     are implied by a zero imaginary part.
-
-      call adjfft(field, nx, ny, nz, in, jn)
-
-      return
-      end
-
-
-c===================================================================
-
-      subroutine processk(i, j, k, dk, PSMin, PSStep, PSTable, 
-     &                    itype, z, kcutoffsq, box)
-
-      implicit none
-
-!     Parameter
-
-      real, parameter :: twopi = 2.0*3.14159265358979324
-
-!     Arguments
-
-      integer :: i, j, k, itype
-      real ::    dk, PSMin, PSStep, PSTable(*)
-      real ::    kcutoffsq, box
-      complex :: z
-
-!     Locals
-
-      integer :: index
-      real :: psval, kdir, klog, ang, amp, kmodsq
-      real :: ranf_min
-
-!     External function
-
-      real :: ran1
-      real :: enzo_ranf
-
-
-
-!     Define table lookup function
- 
-      real ::    Table1, Table2, Step, Min, Tablex, TableLookUp
-      integer :: Tablei
-
-      TableLookUp(Table1, Table2, Step, Min, Tablei, Tablex) = 
-     &            Table1 + (Tablex - real(Tablei-1)*Step - Min) 
-     &            / Step * (Table2 - Table1)
-
-
-      kmodsq = max(i**2 + j**2 + k**2, 1)*dk**2
-      klog   = 0.5*log(kmodsq)
-      index = int((klog - PSMin)/PSStep)
-      psval = TableLookUp(PSTable(index), PSTable(index+1),
-     &                    PSStep, PSMin, index, klog)
-      psval = psval * dk**3
-
-      if (kmodsq .gt. kcutoffsq) psval = 0.0
-
-!     Generate a complex number with random phase and amplitude
-!     Gaussian distributed with a mean of sqrt(psval) with the
-!     Box-Muller method.  Note we have supressed a factor of
-!     sqrt(2) since we must also divide by this factor to account
-!     for the dreary fact that we are really generating two random
-!     fields (if we were doing a complex-to-complex transform
-!     this would show up when we discarded the perfectly
-!     good imaginary component of the transformed field).  whew.
-
-      ranf_min = 1.0e-37
-
-      ang = twopi*enzo_ranf()
-      amp = sqrt(-log(max(enzo_ranf(),ranf_min)) * psval)
-      z   = cmplx(cos(ang), sin(ang)) * amp
-
-!     Process this on the basis of itype:
-!      0)   density field - just leave it be.
-!      1-3) displacement field - multiply by vec(k)/k^2
-!           (and then convert from Mpc to fraction of box).
-
-      if (itype .ne. 0) then
-         if (itype .eq. 1) kdir = real(i)*dk
-         if (itype .eq. 2) kdir = real(j)*dk
-         if (itype .eq. 3) kdir = real(k)*dk
-         z = z * cmplx(0.0,1.0) * kdir / (kmodsq * box)
-      endif
-
-      return
-      end

File src/inits/acml_st1.F

View file
  • Ignore whitespace
+#include "../enzo/fortran.def"
+
+#ifdef XT3
+
+#ifdef CONFIG_BFLOAT_4
+
+      subroutine acml_st1(x, n1, idir)
+
+      implicit none
+#include "../enzo/fortran_types.def"
+
+      INTG_PREC :: n1, idir
+      CMPLX_PREC :: x(n1)
+
+      REAL*4 :: factor
+      REAL*4 :: scale
+      complex*8, allocatable :: work(:)
+
+      integer*4 :: nwork, jdir
+      integer*4 :: m1, info, i
+
+      m1 = n1
+      nwork = 5*n1+100
+      jdir = idir
+
+      allocate( work(nwork) )
+
+      call cfft1d(0_IKIND, m1, x, work, info)
+      call cfft1d(jdir,    m1, x, work, info)
+
+      deallocate( work )
+
+      if( idir == -1 ) then
+        do i = 1, n1
+          x(i) = x(i) * sqrt(REAL(n1,RKIND))
+        end do
+      else
+        do i = 1, n1
+          x(i) = x(i) / sqrt(REAL(n1,RKIND))
+        end do
+      end if
+
+      return
+      end
+
+#endif
+
+#ifdef CONFIG_BFLOAT_8
+
+      subroutine acml_st1(x, n1, idir)
+
+      implicit none
+#include "../enzo/fortran_types.def"
+
+      INTG_PREC :: n1, idir
+      CMPLX_PREC :: x(n1)
+
+      REAL*8 :: factor
+      REAL*8 :: scale
+      complex*16, allocatable :: work(:)
+
+      integer*4 :: nwork, jdir
+      integer*4 :: m1, info, i
+
+      m1 = n1
+      nwork = 5*n1+100
+      jdir = idir
+
+      allocate( work(nwork) )
+
+      call zfft1d(0_IKIND, m1, x, work, info)
+      if( info .ne. 0 ) write(0,'("Info1 = ",i4)') info
+      call zfft1d(jdir,    m1, x, work, info)
+      if( info .ne. 0 ) write(0,'("Info2 = ",i4)') info
+
+
+      deallocate( work )
+
+      if( idir == -1 ) then
+        do i = 1, n1
+          x(i) = x(i) * sqrt(REAL(n1,RKIND))
+        end do
+      else
+        do i = 1, n1
+          x(i) = x(i) / sqrt(REAL(n1,RKIND))
+        end do
+      end if
+
+      return
+      end
+
+#endif
+
+#else
+
+      subroutine acml_st1(x, n1, idir)
+
+      implicit none
+#include "../enzo/fortran_types.def"
+
+      INTG_PREC n1, idir
+      CMPLX_PREC x(n1)
+
+      write(0,'("ACML stride 1 FFT error")')
+      call stop_all_cpus
+
+      return
+      end
+
+#endif

File src/inits/acml_st1.src

  • Ignore whitespace
-#ifdef XT3
-
-#ifdef CONFIG_BFLOAT_4
-
-      subroutine acml_st1(x, n1, idir)
-
-      implicit none
-
-      integer :: n1, idir
-      complex :: x(n1)
-
-      real*4 :: factor
-      real*4 :: scale
-      complex*8, allocatable :: work(:)
-
-      integer*4 :: nwork, jdir
-      integer*4 :: m1, info, i
-
-      m1 = n1
-      nwork = 5*n1+100
-      jdir = idir
-
-      allocate( work(nwork) )
-
-      call cfft1d(   0, m1, x, work, info)
-      call cfft1d(jdir, m1, x, work, info)
-
-      deallocate( work )
-
-      if( idir == -1 ) then
-        do i = 1, n1
-          x(i) = x(i) * sqrt(real(n1))
-        end do
-      else
-        do i = 1, n1
-          x(i) = x(i) / sqrt(real(n1))
-        end do
-      end if
-
-      return
-      end
-
-#endif
-
-#ifdef CONFIG_BFLOAT_8
-
-      subroutine acml_st1(x, n1, idir)
-
-      implicit none
-
-      integer :: n1, idir
-      complex :: x(n1)
-
-      real*8 :: factor
-      real*8 :: scale
-      complex*16, allocatable :: work(:)
-
-      integer*4 :: nwork, jdir
-      integer*4 :: m1, info, i
-
-      m1 = n1
-      nwork = 5*n1+100
-      jdir = idir
-
-      allocate( work(nwork) )
-
-      call zfft1d(   0, m1, x, work, info)
-      if( info .ne. 0 ) write(0,'("Info1 = ",i4)') info
-      call zfft1d(jdir, m1, x, work, info)
-      if( info .ne. 0 ) write(0,'("Info2 = ",i4)') info
-
-
-      deallocate( work )
-
-      if( idir == -1 ) then
-        do i = 1, n1
-          x(i) = x(i) * sqrt(real(n1))
-        end do
-      else
-        do i = 1, n1
-          x(i) = x(i) / sqrt(real(n1))
-        end do
-      end if
-
-      return
-      end
-
-#endif
-
-#else
-
-      subroutine acml_st1(x, n1, idir)
-
-      implicit none
-
-      integer :: n1, idir
-      complex :: x(n1)
-
-      write(0,'("ACML stride 1 FFT error")')
-      call stop_all_cpus
-
-      return
-      end
-
-#endif

File src/inits/adjfft.F

View file
  • Ignore whitespace
+#include "../enzo/fortran.def"
+!=======================================================================
+!////////////////////////  SUBROUTINE ADJFFT  \\\\\\\\\\\\\\\\\\\\\\\\\\
+
+      subroutine adjfft(array, nxz, nyz, nzz, in, jn)
+
+!  ADJUSTS A COMPLEX ARRAY TO SATISFY THE CONJUGATE RELATIONS FOR REAL ARRAYS
+!
+!  written by: Greg Bryan
+!  date:       January, 1993
+!  modified:   Robert Harkness
+!  date:       November, 2003
+!
+!  INPUTS: array(in,jn,*) - 3d array to be adjusted
+!          in,jn        - first and second declared dimension of array
+!          nxz,nyz,nzz  - size (within array) of appropriate 3d array
+!
+!  OUTPUTS: array(in,jn,*)
+
+      implicit NONE
+#include "../enzo/fortran_types.def"
+
+!     Arguments
+
+      INTG_PREC :: nxz, nyz, nzz, in, jn
+      R_PREC :: array(in, jn, *)
+
+!     Locals
+
+      INTG_PREC :: j, k
+
+!-----------------------------------------------------------------------
+
+!  Adjust the array to satisfy the complex conjugation relations for real
+!  arrays.
+
+!  1) adjust corners
+
+      array(    2,       1,       1) = 0.0_RKIND
+      array(nxz+2,       1,       1) = 0.0_RKIND
+      array(    2, nyz/2+1,       1) = 0.0_RKIND
+      array(    2,       1, nzz/2+1) = 0.0_RKIND
+      array(nxz+2, nyz/2+1,       1) = 0.0_RKIND
+      array(nxz+2,       1, nzz/2+1) = 0.0_RKIND
+      array(    2, nyz/2+1, nzz/2+1) = 0.0_RKIND
+      array(nxz+2, nyz/2+1, nzz/2+1) = 0.0_RKIND
+
+!  2) adjust faces
+
+      do j = 2, nyz/2
+         array(    1, j,       1) =  array(    1, nyz+2-j,       1)
+         array(    2, j,       1) = -array(    2, nyz+2-j,       1)
+         array(nxz+1, j,       1) =  array(nxz+1, nyz+2-j,       1)
+         array(nxz+2, j,       1) = -array(nxz+2, nyz+2-j,       1)
+         array(    1, j, nzz/2+1) =  array(    1, nyz+2-j, nzz/2+1)
+         array(    2, j, nzz/2+1) = -array(    2, nyz+2-j, nzz/2+1)
+         array(nxz+1, j, nzz/2+1) =  array(nxz+1, nyz+2-j, nzz/2+1)
+         array(nxz+2, j, nzz/2+1) = -array(nxz+2, nyz+2-j, nzz/2+1)
+         do k = 2, nzz/2
+            array(    1, j,       k) =  array(    1, nyz+2-j, nzz+2-k)
+            array(    2, j,       k) = -array(    2, nyz+2-j, nzz+2-k)
+            array(nxz+1, j,       k) =  array(nxz+1, nyz+2-j, nzz+2-k)
+            array(nxz+2, j,       k) = -array(nxz+2, nyz+2-j, nzz+2-k)
+            array(    1, j, nzz+2-k) =  array(    1, nyz+2-j,       k)
+            array(    2, j, nzz+2-k) = -array(    2, nyz+2-j,       k)
+            array(nxz+1, j, nzz+2-k) =  array(nxz+1, nyz+2-j,       k)
+            array(nxz+2, j, nzz+2-k) = -array(nxz+2, nyz+2-j,       k)
+         enddo
+      enddo
+
+!  3) adjust sides
+
+      do k = 2, nzz/2
+         array(    1,       1, k) =  array(    1,       1, nzz+2-k)
+         array(    2,       1, k) = -array(    2,       1, nzz+2-k)
+         array(nxz+1,       1, k) =  array(nxz+1,       1, nzz+2-k)
+         array(nxz+2,       1, k) = -array(nxz+2,       1, nzz+2-k)
+         array(    1, nyz/2+1, k) =  array(    1, nyz/2+1, nzz+2-k)
+         array(    2, nyz/2+1, k) = -array(    2, nyz/2+1, nzz+2-k)
+         array(nxz+1, nyz/2+1, k) =  array(nxz+1, nyz/2+1, nzz+2-k)
+         array(nxz+2, nyz/2+1, k) = -array(nxz+2, nyz/2+1, nzz+2-k)
+      enddo
+
+      return
+      end

File src/inits/adjfft.src

  • Ignore whitespace
-!=======================================================================
-!////////////////////////  SUBROUTINE ADJFFT  \\\\\\\\\\\\\\\\\\\\\\\\\\
-
-      subroutine adjfft(array, nxz, nyz, nzz, in, jn)
-
-!  ADJUSTS A COMPLEX ARRAY TO SATISFY THE CONJUGATE RELATIONS FOR REAL ARRAYS
-!
-!  written by: Greg Bryan
-!  date:       January, 1993
-!  modified:   Robert Harkness
-!  date:       November, 2003
-!
-!  INPUTS: array(in,jn,*) - 3d array to be adjusted
-!          in,jn        - first and second declared dimension of array
-!          nxz,nyz,nzz  - size (within array) of appropriate 3d array
-!
-!  OUTPUTS: array(in,jn,*)
-
-      use enzo_precision
-
-      implicit NONE
-
-!     Arguments
-
-      integer (enzo_int) :: nxz, nyz, nzz, in, jn
-      real    (enzo_fpr) :: array(in, jn, *)
-
-!     Locals
-
-      integer :: j, k
-
-!-----------------------------------------------------------------------
-
-!  Adjust the array to satisfy the complex conjugation relations for real
-!  arrays.
-
-!  1) adjust corners
-
-      array(    2,       1,       1) = 0.0
-      array(nxz+2,       1,       1) = 0.0
-      array(    2, nyz/2+1,       1) = 0.0
-      array(    2,       1, nzz/2+1) = 0.0
-      array(nxz+2, nyz/2+1,       1) = 0.0
-      array(nxz+2,       1, nzz/2+1) = 0.0
-      array(    2, nyz/2+1, nzz/2+1) = 0.0
-      array(nxz+2, nyz/2+1, nzz/2+1) = 0.0
-
-!  2) adjust faces
-
-      do j = 2, nyz/2
-         array(    1, j,       1) =  array(    1, nyz+2-j,       1)
-         array(    2, j,       1) = -array(    2, nyz+2-j,       1)
-         array(nxz+1, j,       1) =  array(nxz+1, nyz+2-j,       1)
-         array(nxz+2, j,       1) = -array(nxz+2, nyz+2-j,       1)
-         array(    1, j, nzz/2+1) =  array(    1, nyz+2-j, nzz/2+1)
-         array(    2, j, nzz/2+1) = -array(    2, nyz+2-j, nzz/2+1)
-         array(nxz+1, j, nzz/2+1) =  array(nxz+1, nyz+2-j, nzz/2+1)
-         array(nxz+2, j, nzz/2+1) = -array(nxz+2, nyz+2-j, nzz/2+1)
-         do k = 2, nzz/2
-            array(    1, j,       k) =  array(    1, nyz+2-j, nzz+2-k)
-            array(    2, j,       k) = -array(    2, nyz+2-j, nzz+2-k)
-            array(nxz+1, j,       k) =  array(nxz+1, nyz+2-j, nzz+2-k)
-            array(nxz+2, j,       k) = -array(nxz+2, nyz+2-j, nzz+2-k)
-            array(    1, j, nzz+2-k) =  array(    1, nyz+2-j,       k)
-            array(    2, j, nzz+2-k) = -array(    2, nyz+2-j,       k)
-            array(nxz+1, j, nzz+2-k) =  array(nxz+1, nyz+2-j,       k)
-            array(nxz+2, j, nzz+2-k) = -array(nxz+2, nyz+2-j,       k)
-         enddo
-      enddo
-
-!  3) adjust sides
-
-      do k = 2, nzz/2
-         array(    1,       1, k) =  array(    1,       1, nzz+2-k)
-         array(    2,       1, k) = -array(    2,       1, nzz+2-k)
-         array(nxz+1,       1, k) =  array(nxz+1,       1, nzz+2-k)
-         array(nxz+2,       1, k) = -array(nxz+2,       1, nzz+2-k)
-         array(    1, nyz/2+1, k) =  array(    1, nyz/2+1, nzz+2-k)
-         array(    2, nyz/2+1, k) = -array(    2, nyz/2+1, nzz+2-k)
-         array(nxz+1, nyz/2+1, k) =  array(nxz+1, nyz/2+1, nzz+2-k)
-         array(nxz+2, nyz/2+1, k) = -array(nxz+2, nyz/2+1, nzz+2-k)
-      enddo
-
-      return
-      end

File src/inits/cosmo_functions.F

View file
  • Ignore whitespace
+#include "../enzo/fortran.def"
+!=======================================================================
+! Set the current radiation density (incl. massless neutrinos) for h=1
+!define OMEGA0_RAD 4.22e-5
+#define OMEGA0_RAD 0.0
+!=======================================================================
+
+!=======================================================================
+! Set R_PREC to match what's in fortran_types.def (since it cannot be 
+! included before the first function declaration)
+#ifdef CONFIG_BFLOAT_4
+#define R_PREC real*4
+#endif
+
+#ifdef CONFIG_BFLOAT_8
+#define R_PREC real*8
+#endif
+!=======================================================================
+
+
+!=======================================================================
+!//////////////////////  FUNCTION COMPUTE_TIME  \\\\\\\\\\\\\\\\\\\\\\\\
+
+      R_PREC function calc_time(aye_temp)
+
+!  COMPUTES THE TIME GIVEN THE EXPANSION FACTOR
+!
+!  written by: Greg Bryan
+!  date:       February, 1992
+!  modified:   Robert Harkness
+!  date:       November, 2003
+!
+!  PURPOSE:  
+!
+!  INPUTS:
+!
+!  OUTPUTS:
+
+      implicit none
+#include "../enzo/fortran_types.def"
+
+#include "cosmo.h"
+
+!     Argument
+
+      R_PREC :: aye_temp
+
+!     Locals
+
+      R_PREC :: time_temp
+
+!     Externals
+
+      R_PREC :: dtda, midpnt
+      external :: dtda, midpnt
+
+
+      call qromo(dtda, 0.0_RKIND, aye_temp, time_temp, midpnt)
+
+      calc_time = time_temp
+
+      return
+      end
+
+
+!=======================================================================
+!/////////////////////////  FUNCTION DADT  \\\\\\\\\\\\\\\\\\\\\\\\\\\\\
+
+      R_PREC function calc_ayed(aye_temp)
+
+!  modified:   Robert Harkness
+!  date:       November, 2003
+
+      implicit none
+#include "../enzo/fortran_types.def"
+
+!     Argument
+
+      R_PREC :: aye_temp
+
+!     External function
+
+      R_PREC :: dtda
+
+
+      calc_ayed = 1.0_RKIND/dtda(aye_temp)
+
+      return
+      end
+
+
+!=======================================================================
+!////////////////////////  FUNCTION D2A/DT2  \\\\\\\\\\\\\\\\\\\\\\\\\\\
+
+      R_PREC function calc_ayedd(aye_temp)
+
+!  modified:   Robert Harkness
+!  date:       November, 2003
+
+      implicit none
+#include "../enzo/fortran_types.def"
+
+#include "cosmo.h"
+
+!     Argument
+
+      R_PREC :: aye_temp
+
+!     Locals
+
+      R_PREC :: omega0_rad, omega0_mrad
+
+      omega0_rad = OMEGA0_RAD*hub**(-2)
+      omega0_mrad = omega0 - omega0_rad
+
+      calc_ayedd = aye_temp*(lam0 
+     &     - 0.5_RKIND*omega0_mrad/(aye_temp*uaye)**3
+     &     - omega0_rad /(aye_temp*uaye)**4)
+
+!     Convert to code units (the factors of H_0 have already been cancelled)
+
+     &     /(sqrt(1.5_RKIND*omega0)*(1.0_RKIND+zri)**1.5_RKIND)**2
+
+      return
+      end
+
+
+!=======================================================================
+!/////////////////////////  FUNCTION DTDA  \\\\\\\\\\\\\\\\\\\\\\\\\\\\\
+
+      R_PREC function dtda(aye_temp)
+
+!  COMPUTES THE VALUE OF DT/DA, GIVEN A
+
+!  written by: Greg Bryan
+!  date:       February, 1992
+!  modified:   Robert Harkness
+!  date:       November, 2003
+
+      implicit none
+#include "../enzo/fortran_types.def"
+
+#include "cosmo.h"
+
+!     Argument
+
+      R_PREC :: aye_temp
+
+!     Locals
+
+      R_PREC :: at2, omega0_rad, omega0_mrad
+
+
+!     We include the (small) effect of radiation (Peebles 6.80)
+
+      omega0_rad = OMEGA0_RAD*hub**(-2)
+      omega0_mrad = omega0 - omega0_rad
+
+!     Convert aye from code units
+
+      at2 = aye_temp*uaye
+
+!     Compute dt/da (Peebles1993, p. 312)
+
+      dtda = 1._RKIND/sqrt(omega0_mrad/at2 + omega0_rad/at2**2 +
+     &                lam0*at2**2 + 1._RKIND - lam0 - omega0)
+
+!     Convert to code units (the factors of H_0 have already been cancelled)
+
+     &       *sqrt(1.5_RKIND*omega0)*(1._RKIND+zri)**1.5_RKIND * uaye
+
+      return
+      end
+
+
+!=======================================================================
+!////////////////////////  FUNCTION CALC_AYE  \\\\\\\\\\\\\\\\\\\\\\\\\\
+
+      R_PREC function calc_aye(time_temp)
+
+!  COMPUTES THE EXPANSION FACTOR GIVEN THE TIME
+!
+!  written by: Greg Bryan
+!  date:       February, 1992
+!  modified:   Robert Harkness
+!  date:       November, 2003
+
+      implicit none
+#include "../enzo/fortran_types.def"
+
+#include "cosmo.h"
+
+!     Argument
+
+      R_PREC :: time_temp
+
+!     Locals
+
+      R_PREC :: aye_temp, aye_old, calc_time, dtda, tfinal, tfromfinal
+      INTG_PREC :: i
+
+!     Parameters
+
+      INTG_PREC, parameter :: niter = 10
+      R_PREC, parameter :: tolerance = 1.e-5_RKIND
+
+
+!     Make an initial guess based on Taylor expansion (i.e. use q0)
+
+      tfinal = calc_time(1.0_RKIND+zri)
+      tfromfinal = sqrt(2._RKIND/3._RKIND/omega0) * 
+     &     (1._RKIND+zri)**(-1.5_RKIND) * (tfinal - time_temp)
+      aye_temp = (1._RKIND+zri)*(1._RKIND - tfromfinal - 
+     &     0.5_RKIND*(0.5_RKIND*omega0 - lam0)*tfromfinal**2)
+
+!     Do a Newton-Raphson iteration
+
+      do i = 1, niter
+         aye_old = aye_temp
+         aye_temp = aye_old + 1._RKIND/dtda(aye_old) *
+     &                        (time_temp - calc_time(aye_old))
+         if (abs(aye_old-aye_temp)/aye_temp .lt. tolerance) goto 100
+      enddo
+
+      write(0,*) 'NR in calc_aye failed.'
+      stop
+
+!     Done
+
+ 100  continue
+
+      calc_aye = aye_temp
+
+      return
+      end
+
+
+!=======================================================================
+!/////////////////////////  FUNCTION CALC_F  \\\\\\\\\\\\\\\\\\\\\\\\\\\
+
+      R_PREC function calc_f(aye_temp)
+
+!  COMPUTES THE FUNCTION D LOG (DELTA_PLUS) / D LOG (AYE)
+!
+!  written by: Greg Bryan
+!  date:       February, 1995
+!  modified:   Robert Harkness
+!  date:       November, 2003
+
+      implicit none
+#include "../enzo/fortran_types.def"
+
+#include "cosmo.h"
+
+!     Argument
+
+      R_PREC :: aye_temp
+
+!     Locals
+
+      R_PREC :: ayed, ayedd, sum, at2
+
+      R_PREC :: fhelper, calc_ayed, calc_ayedd, midpnt
+
+      external :: fhelper, midpnt
+
+
+!     We calculate f(z) through PEEBLES93, eq 13.81:
+!
+!     f = ayedd*aye/ayed^2 - 1 + 1 / (a^2 E^3 G)
+!
+!        where G = int(z to infifinty) (1+z)/E(z)^3 dz
+!                = int(0 to a) da/(E(a) * a)^3
+!              E = da/dt / H_0 
+
+!     Compute G (using the usual a convention)
+
+      at2 = aye_temp*uaye
+
+      call qromo(fhelper, 0.0_RKIND, at2, sum, midpnt)
+
+      ayed = calc_ayed(aye_temp)
+      ayedd = calc_ayedd(aye_temp)
+
+!     Note that f is dimensionless, specifically all the unusual aye
+!     convention (a=1 at z=zri) cancels for the first term.  The third
+!     term is computed entirely in the usual convention.
+
+      calc_f = ayedd*aye_temp/ayed**2 - 1.0_RKIND + at2*fhelper(at2)/sum
+
+      return
+      end
+
+
+!=======================================================================
+!/////////////////////////  FUNCTION FFUNC  \\\\\\\\\\\\\\\\\\\\\\\\\\\\\
+
+      R_PREC function fhelper(at2)
+
+!  modified:   Robert Harkness
+!  date:       November, 2003
+
+      implicit none
+#include "../enzo/fortran_types.def"
+
+#include "cosmo.h"
+
+!     returns the function 1/(a * E)^3  
+!       (where, here only, a is not aye, i.e. it obeys the usual convention
+!        of a = 1 at z = 0)
+
+!     Argument
+
+      R_PREC :: at2
+
+!     Locals
+
+      R_PREC :: E, omega0_rad, omega0_mrad
+
+
+      omega0_rad = OMEGA0_RAD*hub**(-2)
+      omega0_mrad = omega0 - omega0_rad
+
+      E = sqrt(omega0_mrad/at2**3 + (1.0_RKIND - omega0 - lam0)/at2**2 + 
+     &         lam0 + omega0_rad/at2**4)
+
+      fhelper = 1.0_RKIND/(at2*E)**3
+
+      return
+      end
+
+
+!=======================================================================
+!///////////////////////  FUNCTION CALC_GROWTH  \\\\\\\\\\\\\\\\\\\\\\\\
+
+      R_PREC function calc_growth(z1)
+
+!  COMPUTES THE GROWTH FUNCTION D(z), NORMALIZED AT HIGH REDSHIFT
+!
+!  written by: Greg Bryan
+!  date:       February, 1995
+!  modified:   Robert Harkness
+!  date:       November, 2003
+
+!  PURPOSE:
+!    Note: this function is _not_ normalized to D(z=0) = 1
+!
+!  INPUTS:
+!    z      - redshift function is to be evaluated at
+!    omega0 - matter density ratio at z=0
+!    omega_lam - \Lambda/(3*H_0^2) at z=0
+!
+!  OUTPUTS:
+!    calc_growth - linear growth, normalized to 1/(1+z) at z=infinity
+
+
+      implicit none
+#include "../enzo/fortran_types.def"
+
+#include "cosmo.h"
+
+!     Argument
+
+      R_PREC :: z1
+
+!     Locals
+
+      R_PREC :: a, sum, E
+
+!     Externals
+
+      R_PREC :: fhelper, midpnt
+      external :: fhelper, midpnt
+
+
+!     We calculate D(z) through PEEBLES93, eq 13.78:
+!
+!     D(z) = E(z) * G(z)
+!
+!        where G = 5*omega0/2 * int(z to infinity) (1+z)/E(z)^3 dz
+!                = 5*omega0/2 * int(0 to a) da/(E(a) * a)^3
+!              E = da/dt / H_0
+
+!     Compute G, missing the prefactor
+
+      a = 1.0_RKIND/(1.0_RKIND+z1)
+
+      call qromo(fhelper, 1.e-10_RKIND, a, sum, midpnt)
+
+      E = sqrt(omega0/a**3 + (1._RKIND - omega0 - lam0)/a**2 + lam0)
+
+      calc_growth = 5.0_RKIND*omega0/2.0_RKIND * sum * E
+
+      return
+      end
+
+
+!=======================================================================
+!/////////////////////  SUBROUTINE SET_COMMON  \\\\\\\\\\\\\\\\\\\\\\\\\
+
+      subroutine set_common(lam0_in, omega0_in, zri_in, hub_in)
+
+!  modified:   Robert Harkness
+!  date:       November, 2003
+
+      implicit none
+#include "../enzo/fortran_types.def"
+
+#include "cosmo.h"
+
+!     Arguments
+
+      R_PREC :: lam0_in, omega0_in, zri_in, hub_in
+
+
+      lam0   = lam0_in
+      omega0 = omega0_in
+      zri    = zri_in
+      hub    = hub_in
+      uaye   = 1.0_RKIND/(1.0_RKIND + zri)
+
+      return
+      end

File src/inits/cosmo_functions.src

  • Ignore whitespace
-!=======================================================================
-! Set the current radiation density (incl. massless neutrinos) for h=1
-!define OMEGA0_RAD 4.22e-5
-#define OMEGA0_RAD 0.0
-!=======================================================================
-
-
-!=======================================================================
-!//////////////////////  FUNCTION COMPUTE_TIME  \\\\\\\\\\\\\\\\\\\\\\\\
-
-      real function calc_time(aye_temp)
-
-!  COMPUTES THE TIME GIVEN THE EXPANSION FACTOR
-!
-!  written by: Greg Bryan
-!  date:       February, 1992
-!  modified:   Robert Harkness
-!  date:       November, 2003
-!
-!  PURPOSE:  
-!
-!  INPUTS:
-!
-!  OUTPUTS:
-
-      implicit none
-
-#include "cosmo.h"
-
-!     Argument
-
-      real :: aye_temp
-
-!     Locals
-
-      real :: time_temp
-
-!     Externals
-
-      real :: dtda, midpnt
-      external :: dtda, midpnt
-
-
-      call qromo(dtda, 0.0, aye_temp, time_temp, midpnt)
-
-      calc_time = time_temp
-
-      return
-      end
-
-
-!=======================================================================
-!/////////////////////////  FUNCTION DADT  \\\\\\\\\\\\\\\\\\\\\\\\\\\\\
-
-      real function calc_ayed(aye_temp)
-
-!  modified:   Robert Harkness
-!  date:       November, 2003
-
-      implicit none
-
-!     Argument
-
-      real :: aye_temp
-
-!     External function
-
-      real :: dtda
-
-
-      calc_ayed = 1.0/dtda(aye_temp)
-
-      return
-      end
-
-
-!=======================================================================
-!////////////////////////  FUNCTION D2A/DT2  \\\\\\\\\\\\\\\\\\\\\\\\\\\
-
-      real function calc_ayedd(aye_temp)
-
-!  modified:   Robert Harkness
-!  date:       November, 2003
-
-      implicit none
-
-#include "cosmo.h"
-
-!     Argument
-
-      real :: aye_temp
-
-!     Locals
-
-      real :: omega0_rad, omega0_mrad
-
-      omega0_rad = OMEGA0_RAD*hub**(-2)
-      omega0_mrad = omega0 - omega0_rad
-
-      calc_ayedd = aye_temp*(lam0 - 0.5*omega0_mrad/(aye_temp*uaye)**3
-     &                            -     omega0_rad /(aye_temp*uaye)**4)
-
-!     Convert to code units (the factors of H_0 have already been cancelled)
-
-     &       /(sqrt(3.0/2.0*omega0)*(1.0+zri)**1.5)**2
-
-      return
-      end
-
-
-!=======================================================================
-!/////////////////////////  FUNCTION DTDA  \\\\\\\\\\\\\\\\\\\\\\\\\\\\\
-
-      real function dtda(aye_temp)
-
-!  COMPUTES THE VALUE OF DT/DA, GIVEN A
-
-!  written by: Greg Bryan
-!  date:       February, 1992
-!  modified:   Robert Harkness
-!  date:       November, 2003
-
-      implicit none
-
-#include "cosmo.h"
-
-!     Argument
-
-      real :: aye_temp
-
-!     Locals
-
-      real :: at2, omega0_rad, omega0_mrad
-
-
-!     We include the (small) effect of radiation (Peebles 6.80)
-
-      omega0_rad = OMEGA0_RAD*hub**(-2)
-      omega0_mrad = omega0 - omega0_rad
-
-!     Convert aye from code units
-
-      at2 = aye_temp*uaye
-
-!     Compute dt/da (Peebles1993, p. 312)
-
-      dtda = 1.0/sqrt(omega0_mrad/at2 + omega0_rad/at2**2 +
-     &                lam0*at2**2 + 1.0 - lam0 - omega0)
-
-!     Convert to code units (the factors of H_0 have already been cancelled)
-
-     &       *sqrt(3.0/2.0*omega0)*(1.0+zri)**1.5 * uaye
-
-      return
-      end
-
-
-!=======================================================================
-!////////////////////////  FUNCTION CALC_AYE  \\\\\\\\\\\\\\\\\\\\\\\\\\
-
-      real function calc_aye(time_temp)
-
-!  COMPUTES THE EXPANSION FACTOR GIVEN THE TIME
-!
-!  written by: Greg Bryan
-!  date:       February, 1992
-!  modified:   Robert Harkness
-!  date:       November, 2003
-
-      implicit none
-
-#include "cosmo.h"
-
-!     Argument
-
-      real :: time_temp
-
-!     Locals
-
-      real :: aye_temp, aye_old, calc_time, dtda, tfinal, tfromfinal
-      integer :: i
-
-!     Parameters
-
-      integer, parameter :: niter = 10
-      real, parameter :: tolerance = 1.0e-5
-
-
-!     Make an initial guess based on Taylor expansion (i.e. use q0)
-
-      tfinal = calc_time(1.0+zri)
-      tfromfinal = sqrt(2.0/3.0/omega0)*(1.0+zri)**(-1.5) * 
-     &             (tfinal - time_temp)
-      aye_temp = (1.0+zri)*(1 - tfromfinal - 
-     &                 0.5*(0.5*omega0 - lam0)*tfromfinal**2)
-
-!     Do a Newton-Raphson iteration
-
-      do i = 1, niter
-         aye_old = aye_temp
-         aye_temp = aye_old + 1.0/dtda(aye_old) *
-     &                        (time_temp - calc_time(aye_old))
-         if (abs(aye_old-aye_temp)/aye_temp .lt. tolerance) goto 100
-      enddo
-
-      write(0,*) 'NR in calc_aye failed.'
-      stop
-
-!     Done
-
- 100  continue
-
-      calc_aye = aye_temp
-
-      return
-      end
-
-
-!=======================================================================
-!/////////////////////////  FUNCTION CALC_F  \\\\\\\\\\\\\\\\\\\\\\\\\\\
-
-      real function calc_f(aye_temp)
-
-!  COMPUTES THE FUNCTION D LOG (DELTA_PLUS) / D LOG (AYE)
-!
-!  written by: Greg Bryan
-!  date:       February, 1995
-!  modified:   Robert Harkness
-!  date:       November, 2003
-
-      implicit none
-
-#include "cosmo.h"
-
-!     Argument
-
-      real :: aye_temp
-
-!     Locals
-
-      real :: ayed, ayedd, sum, at2
-
-      real :: fhelper, calc_ayed, calc_ayedd, midpnt
-
-      external :: fhelper, midpnt
-
-