Commits

Daniel Reynolds  committed 75dd897

Updates to fortran precision handling. I've narrowed the problem to the integer precision handling (real numbers seem to be doing fine), but I'm a bit stuck on how to proceed.

  • Participants
  • Parent commits 3179d54

Comments (0)

Files changed (42)

File src/enzo/BlockSolve.F

          INFO = 4
       ELSE IF( K  .LT.0               )THEN
          INFO = 5
-      ELSE IF( LDA.LT.MAX( 1_IKIND, NROWA ) )THEN
+      ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
          INFO = 8
-      ELSE IF( LDB.LT.MAX( 1_IKIND, NROWB ) )THEN
+      ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN
          INFO = 10
-      ELSE IF( LDC.LT.MAX( 1_IKIND, M     ) )THEN
+      ELSE IF( LDC.LT.MAX( 1, M     ) )THEN
          INFO = 13
       END IF
       IF( INFO.NE.0 )THEN
          INFO = 2
       ELSE IF( N.LT.0 )THEN
          INFO = 3
-      ELSE IF( LDA.LT.MAX( 1_IKIND, M ) )THEN
+      ELSE IF( LDA.LT.MAX( 1, M ) )THEN
          INFO = 6
       ELSE IF( INCX.EQ.0 )THEN
          INFO = 8
          INFO = 5
       ELSE IF( INCY.EQ.0 )THEN
          INFO = 7
-      ELSE IF( LDA.LT.MAX( 1_IKIND, M ) )THEN
+      ELSE IF( LDA.LT.MAX( 1, M ) )THEN
          INFO = 9
       END IF
       IF( INFO.NE.0 )THEN
          INFO = -1
       ELSE IF( N.LT.0 ) THEN
          INFO = -2
-      ELSE IF( LDA.LT.MAX( 1_IKIND, M ) ) THEN
+      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
          INFO = -4
       END IF
       IF( INFO.NE.0 ) THEN
          INFO = -1
       ELSE IF( N.LT.0 ) THEN
          INFO = -2
-      ELSE IF( LDA.LT.MAX( 1_IKIND, M ) ) THEN
+      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
          INFO = -4
       END IF
       IF( INFO.NE.0 ) THEN
          INFO = -2
       ELSE IF( NRHS.LT.0 ) THEN
          INFO = -3
-      ELSE IF( LDA.LT.MAX( 1_IKIND, N ) ) THEN
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
          INFO = -5
-      ELSE IF( LDB.LT.MAX( 1_IKIND, N ) ) THEN
+      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
          INFO = -8
       END IF
       IF( INFO.NE.0 ) THEN
 c
 c       clean-up loop
 c
-   20 m = mod(n,3_IKIND)
+   20 m = mod(n,3)
       if( m .eq. 0 ) go to 40
       do 30 i = 1,m
         dtemp = dx(i)
 c
 c        clean-up loop
 c
-   20 m = mod(n,5_IKIND)
+   20 m = mod(n,5)
       if( m .eq. 0 ) go to 40
       do 30 i = 1,m
         dx(i) = da*dx(i)
          INFO = 5
       ELSE IF( N  .LT.0               )THEN
          INFO = 6
-      ELSE IF( LDA.LT.MAX( 1_IKIND, NROWA ) )THEN
+      ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
          INFO = 9
-      ELSE IF( LDB.LT.MAX( 1_IKIND, M     ) )THEN
+      ELSE IF( LDB.LT.MAX( 1, M     ) )THEN
          INFO = 11
       END IF
       IF( INFO.NE.0 )THEN

File src/enzo/FSProb_RadiationSource.F90

 
         ! get 4 random numbers for each source (3 location, 1 strength)
         call random_number(rnums)
-        i = max(min(int(rnums(4)*Nx,IKIND), Nx-1), 2_IKIND)
-        j = max(min(int(rnums(6)*Ny,IKIND), Ny-1), 2_IKIND)
-        k = max(min(int(rnums(8)*Nz,IKIND), Nz-1), 2_IKIND)
+        i = max(min(int(rnums(4)*Nx,IKIND), Nx-1), 2)
+        j = max(min(int(rnums(6)*Ny,IKIND), Ny-1), 2)
+        k = max(min(int(rnums(8)*Nz,IKIND), Nz-1), 2)
 !        eta(i,j,k) = rnums(10)*h_nu0*REAL(NGammaDot/dV,RKIND)
         eta(i,j,k) = rnums(10)*h_nu0*NGammaDot/dV
 !        print '(A,3(i2,1x),A,es9.2)', '   setting source at ',i,j,k,' with strength ',eta(i,j,k)

File src/enzo/MTLPARAM.h

       R_PREC TEMMIN, DELT, DENMIN, DELD, FREQDEL, FREQMIN
 
       PARAMETER(NUME=12,MAXLN=220)
-      PARAMETER(NIT=200,TEMMIN=3.0,DELT=0.03)
-      PARAMETER(NID=300,DENMIN=-12.0,DELD=0.05,NID2=240)
-      PARAMETER(NIB=400,FREQDEL=0.02,FREQMIN=1.0)
+      PARAMETER(NIT=200,TEMMIN=3._RKIND,DELT=0.03_RKIND)
+      PARAMETER(NID=300,DENMIN=-12._RKIND,DELD=0.05_RKIND,NID2=240)
+      PARAMETER(NIB=400,FREQDEL=0.02_RKIND,FREQMIN=1._RKIND)
 C
       INTG_PREC NJ(12), LLJ(12,30), IPHOT
       R_PREC ABUNJ(12), WJ(12,MAXLN), E3J(12,MAXLN), 

File src/enzo/calc_photo_rates.F

       R_PREC    FNUDEL, AVGNU, FNUMNUH0, FNUMNUHE0, FNUMNUHE20,
      &        QNICK, CROSSS, FJXRY11, FJXRY21,
      &        FJXRY12, FJXRY22, FJXRY13, FJXRY23
-      double precision tbase1, xbase1, dbase1, coolunit, mh
+      real*8 tbase1, xbase1, dbase1, coolunit, mh
 c
 c\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\/////////////////////////////////
 c=======================================================================

File src/enzo/calc_rates.F

       parameter (ireco = 1, ibrco = 1, icico = 1, iceco = 1, 
      &           iphoto_mol = 0)
 
-      double precision everg, evhz, tevk, mh, pi, tiny_cgs, dhuge, kb
+      real*8 everg, evhz, tevk, mh, pi, tiny_cgs, dhuge, kb
       parameter (everg = ev2erg, evhz  = 2.41838d14, 
      &           tevk = 1.1605d4, mh = mass_h, 
      &           pi=pi_val, tiny_cgs = 1.0d-37, 
 c
 c         Set various cutoff values in eV.
 c
-      double precision e24,e25,e26,e27,e28a,e28b,e28c,e29a,e29b,e29c,
+      real*8 e24,e25,e26,e27,e28a,e28b,e28c,e29a,e29b,e29c,
      &                 e30a,e30b
         PARAMETER(
      &    e24  = 13.6d0
 c  Locals
 c
       INTG_PREC i,j,n
-      double precision nu, delnu, hnu, logttt, ttt, tev, logtev, 
+      real*8 nu, delnu, hnu, logttt, ttt, tev, logtev, 
      &                 xx, dum, tbase1, xbase1, kunit, coolunit,
      &                 dbase1, dlogtem, kunit_3bdy, x, y, cierate,
      &                 grain_coef, d_ttt, d_logttt, d_dlogtem,
      &                 ttt2, d_ttt2
-      double precision sigma24(nnu), sigma25(nnu), sigma26(nnu), 
+      real*8 sigma24(nnu), sigma25(nnu), sigma26(nnu), 
      &                 sigma27(nnu), sigma28(nnu), sigma29(nnu), 
      &                 sigma30(nnu), sigma31(nnu), fext(nnu)
-      double precision tm
-      double precision HDLR, HDLV, lt, t3, lt3
+      real*8 tm
+      real*8 HDLR, HDLV, lt, t3, lt3
 c
 c\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\/////////////////////////////////////
 c=======================================================================

File src/enzo/calc_tdust_1d.F

 
       INTG_PREC i, iter, c_done, c_total, nm_done
 
-      double precision pert_i, trad4
+      real*8 pert_i, trad4
 
 !  Slice Locals
 
-      double precision kgr(in), kgrplus(in), sol(in), solplus(in), 
+      real*8 kgr(in), kgrplus(in), sol(in), solplus(in), 
      &     slope(in), tdplus(in), tdustnow(in), tdustold(in), pert(in),
      &     bi_t_mid(in), bi_t_high(in)
       logical nm_itmask(in), bi_itmask(in)
 
-      double precision sola, solb, solpa, solpb
+      real*8 sola, solb, solpa, solpb
 
 !\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\/////////////////////////////////
 !=======================================================================
          do i = is+1, ie+1
             if ( nm_itmask(i) ) then
 
-               tdplus(i) = max(0.001_RKIND, ((1+pert(i)) * tdustnow(i)))
+               tdplus(i) = max(1.e-3_RKIND, ((1._RKIND+pert(i)) 
+     &              * tdustnow(i)))
 
             endif
          enddo
 
       INTG_PREC in, is, ie
       R_PREC t_subl
-      double precision tdust(in)
+      real*8 tdust(in)
 
 !  Iteration mask
 
 
 !  Locals
 
-      integer i
+      INTG_PREC i
 
 !  Slice Locals
 
-      double precision kgr(in)
+      real*8 kgr(in)
 
 !\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\/////////////////////////////////
 !=======================================================================
                kgr(i) = kgr200
             else
                kgr(i) = max(tiny, 
-     &              (kgr200 * (tdust(i) / 1500._RKIND)**-12))
+     &              (kgr200 * (tdust(i) / 1.5e3_RKIND)**(-12)))
             endif
 
          endif

File src/enzo/calc_tdust_3d.F

 
 !  Parameters
 
-      double precision mh
+      real*8 mh
       parameter (mh = mass_h)      !DPC
 
 !  Locals
 
       INTG_PREC i, j, k
       R_PREC trad, zr, logtem0, logtem9, dlogtem
-      double precision coolunit, dbase1, tbase1, xbase1
+      real*8 coolunit, dbase1, tbase1, xbase1
 
 !  Slice locals
  
 !     Compute index into the table and precompute parts of linear interp
 
                indixe(i) = min(nratec-1,
-     &              max(1,int((logtem(i)-logtem0)/dlogtem)+1))
+     &              max(1,int((logtem(i)-logtem0)/dlogtem,IKIND)+1))
                t1(i) = (logtem0 + (indixe(i) - 1)*dlogtem)
                t2(i) = (logtem0 + (indixe(i)    )*dlogtem)
                tdef(i) = (logtem(i) - t1(i)) / (t2(i) - t1(i))

File src/enzo/cic_flag.F

                  ffield(i1  ,1,1) = 1
                  ffield(i1+1,1,1) = 1
               else
-                 istart = max(i1-buffersize+1, 1_IKIND)
+                 istart = max(i1-buffersize+1, 1)
                  iend = min(i1+buffersize, dim1)
                  do ii = istart, iend
                     ffield(ii,1,1) = 1
                  ffield(i1  ,j1+1,1) = 1
                  ffield(i1+1,j1+1,1) = 1
               else
-                 istart = max(i1-buffersize+1, 1_IKIND)
+                 istart = max(i1-buffersize+1, 1)
                  iend = min(i1+buffersize, dim1)
-                 jstart = max(j1-buffersize+1, 1_IKIND)
+                 jstart = max(j1-buffersize+1, 1)
                  jend = min(j1+buffersize, dim2)
                  do jj = jstart, jend
                     do ii = istart, iend
                  ffield(i1  ,j1+1,k1+1) = 1
                  ffield(i1+1,j1+1,k1+1) = 1
               else
-                 istart = max(i1-buffersize+1, 1_IKIND)
+                 istart = max(i1-buffersize+1, 1)
                  iend = min(i1+buffersize, dim1)
-                 jstart = max(j1-buffersize+1, 1_IKIND)
+                 jstart = max(j1-buffersize+1, 1)
                  jend = min(j1+buffersize, dim2)
-                 kstart = max(k1-buffersize+1, 1_IKIND)
+                 kstart = max(k1-buffersize+1, 1)
                  kend = min(k1+buffersize, dim3)
 c	         print*,'cic_flag: buffersize is:',buffersize
                  do kk = kstart, kend

File src/enzo/cluster_maker.F

       R_PREC   pi, G, sndspdC, m1
       R_PREC   isosndsp2, bmass, jeanmass, lum, min_temp
       P_PREC this_pos(3)
-      double precision msolar, mh
+      real*8 msolar, mh
       parameter (pi=pi_val, G=GravConst, 
      &     sndspdC=1.3095e8_RKIND, msolar=SolarMass, mh=mass_h)
 c

File src/enzo/colh2diss.F

 c  Arguments
 c
       R_PREC f1,f2,f3,f4,f5,f6,f7
-      double precision t, tgas
+      real*8 t, tgas
 c
 c  Parameters
 c
 c
 c  Locals
 c
-      double precision tl, a1, a, b, b1, c, c1, d, ftd
-      double precision y(21)
+      real*8 tl, a1, a, b, b1, c, c1, d, ftd
+      real*8 y(21)
 c
 c\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\/////////////////////////////////
 c=======================================================================

File src/enzo/coll_rates.F

 c
       R_PREC k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, k12, k13,
      $       k14, k15, k16, k17, k18, k19, k23
-      double precision kunit, T
+      real*8 kunit, T
       INTG_PREC casebrates
 c
 c  Parameters
 c  Locals
 c
       INTG_PREC i
-      double precision log_T, log_T_eV, T_ev, xx, dum
+      real*8 log_T, log_T_eV, T_ev, xx, dum
 c
 c\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\/////////////////////////////////
 c=======================================================================

File src/enzo/cool1d.F

 c
       R_PREC pmin
       parameter (pmin = tiny)
-      double precision mh
+      real*8 mh
 
       parameter (mh = mass_h)  ! DPC
 
 c
 c        Compute index into the table and precompute parts of linear interp
 c
-         indixe(i) = min(nratec-1, max(1_IKIND,
-     &        int((logtem(i)-logtem0)/dlogtem,IKIND)+1_IKIND))
+         indixe(i) = min(nratec-1, max(1,
+     &        int((logtem(i)-logtem0)/dlogtem,IKIND)+1))
          t1(i) = (logtem0 + (indixe(i) - 1)*dlogtem)
          t2(i) = (logtem0 + (indixe(i)    )*dlogtem)
          tdef(i) = t2(i) - t1(i)

File src/enzo/cool1d_cloudy.F

       R_PREC    comp2, dom, zr
       R_PREC    d(in,jn,kn), de(in,jn,kn), rhoH(in), metallicity(in), 
      &     logtem(in)
-      double precision edot(in)
+      real*8 edot(in)
 
 !  Cloudy parameters and data
 
 
 !     Calculate interpolation index
 
-      index1 = min(gridDim(1)-1, max(1_IKIND,
-     &     int((input1-gridPar1(1))/dgridPar1,IKIND)+1_IKIND))
+      index1 = min(gridDim(1)-1, max(1,
+     &     int((input1-gridPar1(1))/dgridPar1,IKIND)+1))
 
 !     Interpolate over parameter 1
 
 
 !     Calculate interpolation indices
 
-      index1 = min(gridDim(1)-1, max(1_IKIND,
-     &     int((input1-gridPar1(1))/dgridPar1,IKIND)+1_IKIND))
-      index2 = min(gridDim(2)-1, max(1_IKIND,
-     &     int((input2-gridPar2(1))/dgridPar2,IKIND)+1_IKIND))
+      index1 = min(gridDim(1)-1, max(1,
+     &     int((input1-gridPar1(1))/dgridPar1,IKIND)+1))
+      index2 = min(gridDim(2)-1, max(1,
+     &     int((input2-gridPar2(1))/dgridPar2,IKIND)+1))
 
       do q=1, 2
 
 
 !     Calculate interpolation indices
 
-      index1 = min(gridDim(1)-1_IKIND, max(1_IKIND,
-     &     int((input1-gridPar1(1))/dgridPar1,IKIND)+1_IKIND))
-      index2 = min(gridDim(2)-1_IKIND, max(1_IKIND,
-     &     int((input2-gridPar2(1))/dgridPar2,IKIND)+1_IKIND))
-      index3 = min(gridDim(3)-1_IKIND, max(1_IKIND,
-     &     int((input3-gridPar3(1))/dgridPar3,IKIND)+1_IKIND))
+      index1 = min(gridDim(1)-1, max(1,
+     &     int((input1-gridPar1(1))/dgridPar1,IKIND)+1))
+      index2 = min(gridDim(2)-1, max(1,
+     &     int((input2-gridPar2(1))/dgridPar2,IKIND)+1))
+      index3 = min(gridDim(3)-1, max(1,
+     &     int((input3-gridPar3(1))/dgridPar3,IKIND)+1))
 
       do q=1, 2
 
 
 !     Calculate interpolation indices
 
-      index1 = min(gridDim(1)-1_IKIND, max(1_IKIND,
-     &     int((input1-gridPar1(1))/dgridPar1,IKIND)+1_IKIND))
-      index2 = min(gridDim(2)-1_IKIND, max(1_IKIND,
-     &     int((input2-gridPar2(1))/dgridPar2,IKIND)+1_IKIND))
-      index3 = min(gridDim(3)-1_IKIND, max(1_IKIND,
-     &     int((input3-gridPar3(1))/dgridPar3,IKIND)+1_IKIND))
-      index4 = min(gridDim(4)-1_IKIND, max(1_IKIND,
-     &     int((input4-gridPar4(1))/dgridPar4,IKIND)+1_IKIND))
+      index1 = min(gridDim(1)-1, max(1,
+     &     int((input1-gridPar1(1))/dgridPar1,IKIND)+1))
+      index2 = min(gridDim(2)-1, max(1,
+     &     int((input2-gridPar2(1))/dgridPar2,IKIND)+1))
+      index3 = min(gridDim(3)-1, max(1,
+     &     int((input3-gridPar3(1))/dgridPar3,IKIND)+1))
+      index4 = min(gridDim(4)-1, max(1,
+     &     int((input4-gridPar4(1))/dgridPar4,IKIND)+1))
 
       do q=1, 2
 
 
 !     Calculate interpolation indices
 
-      index1 = min(gridDim(1)-1_IKIND, max(1_IKIND,
-     &     int((input1-gridPar1(1))/dgridPar1,IKIND)+1_IKIND))
-      index2 = min(gridDim(2)-1_IKIND, max(1_IKIND,
-     &     int((input2-gridPar2(1))/dgridPar2,IKIND)+1_IKIND))
-      index3 = min(gridDim(3)-1_IKIND, max(1_IKIND,
-     &     int((input3-gridPar3(1))/dgridPar3,IKIND)+1_IKIND))
+      index1 = min(gridDim(1)-1, max(1,
+     &     int((input1-gridPar1(1))/dgridPar1,IKIND)+1))
+      index2 = min(gridDim(2)-1, max(1,
+     &     int((input2-gridPar2(1))/dgridPar2,IKIND)+1))
+      index3 = min(gridDim(3)-1, max(1,
+     &     int((input3-gridPar3(1))/dgridPar3,IKIND)+1))
 #define INDEX_4_BISECTION
 #ifdef INDEX_4_BISECTION
 !     get index 4 with bisection, since not evenly spaced
          enddo
       endif
 #else
-      index4 = min(gridDim(4)-1_IKIND, max(1_IKIND,
-     &     int((input4-gridPar4(1))/dgridPar4,IKIND)+1_IKIND))
+      index4 = min(gridDim(4)-1, max(1,
+     &     int((input4-gridPar4(1))/dgridPar4,IKIND)+1))
 #endif /* INDEX_4_BISECTION */
-      index5 = min(gridDim(5)-1_IKIND, max(1_IKIND,
-     &     int((input5-gridPar5(1))/dgridPar5,IKIND)+1_IKIND))
+      index5 = min(gridDim(5)-1, max(1,
+     &     int((input5-gridPar5(1))/dgridPar5,IKIND)+1))
 
       do q=1, 2
 

File src/enzo/cool1d_koyama.F

 c
       R_PREC pmin
       parameter (pmin = tiny)
-      double precision mh
+      real*8 mh
 
       parameter (mh = mass_h)  ! DPC
 

File src/enzo/cool1d_multi.F

 
 !  Parameters
 
-      double precision mh
+      real*8 mh
       R_PREC    ZSOLAR, mu_metal, boost_factor
 ! MKRJ 12/21/07 boost gammaha by this factor (JHW reduced by 1.0 because
 ! it has been included in v2.0
       R_PREC dt2, ttmin, comp1, comp2, scoef, acoef, HIdot,
      &     hdlte1, hdlow1, gamma2, x, fudge, fH2,
      &     gphdl1, factor, dom_inv, tau, ciefudge
-      double precision coolunit, dbase1, tbase1, xbase1, rtunits,
+      real*8 coolunit, dbase1, tbase1, xbase1, rtunits,
      &     nH2, nother
 
 !  Slice locals
       R_PREC t1(in), t2(in), logtem(in), tdef(in), p2d(in),
      &     tgas(in), tgasold(in), tdust(in), rhoH(in), nh(in),
      &     metallicity(in)
-      double precision edot(in)
+      real*8 edot(in)
 
 !  Cooling/heating slice locals
 
-      double precision ceHI(in), ceHeI(in), ceHeII(in),
+      real*8 ceHI(in), ceHeI(in), ceHeII(in),
      &     ciHI(in), ciHeI(in), ciHeIS(in), ciHeII(in),
      &     reHII(in), reHeII1(in), reHeII2(in), reHeIII(in),
      &     brem(in), cieco(in)
 
       logtem0 = log(temstart)
       logtem9 = log(temend)
-      dlogtem= (log(temend) - log(temstart))/REAL(nratec-1,RKIND)
+      dlogtem= (log(temend) - log(temstart))/REAL(nratec-1)
 
 !     Set units
 
                p2d(i) = e(i,j,k) - 0.5_RKIND*u(i,j,k)**2
                if (idim .gt. 1) p2d(i) = p2d(i) - 0.5_RKIND*v(i,j,k)**2
                if (idim .gt. 2) p2d(i) = p2d(i) - 0.5_RKIND*w(i,j,k)**2
-               p2d(i) = max((gamma - 1._RKIND)*d(i,j,k)*p2d(i), 
-     &              REAL(tiny,RKIND))
+               p2d(i) = max((gamma - 1._RKIND)*d(i,j,k)*p2d(i), tiny)
                end if
             enddo
          endif
 
 !        Compute index into the table and precompute parts of linear interp
 
-         indixe(i) = min(nratec-1, max(1_IKIND,
-     &        int((logtem(i)-logtem0)/dlogtem,IKIND)+1_IKIND))
+         indixe(i) = min(nratec-1, max(1,
+     &        int((logtem(i)-logtem0)/dlogtem,IKIND)+1))
          t1(i) = (logtem0 + (indixe(i) - 1)*dlogtem)
          t2(i) = (logtem0 + (indixe(i)    )*dlogtem)
          tdef(i) = (logtem(i) - t1(i)) / (t2(i) - t1(i))
 #endif /* OPTICAL_DEPTH_FUDGE */
 
             edot(i) = edot(i) - float(ih2co)*fudge*H2I(i,j,k)*(
-     &           vibh(i)/(1._RKIND+vibh(i)/max(   vibl     ,tiny)) +
+     &           vibh(i)/(1._RKIND+vibh(i)/max(   vibl   ,tiny)) +
      &           roth(i)/(1._RKIND+roth(i)/max(qq*rotl(i),tiny))     
      &           )/2._RKIND/dom
             end if
             do i = is+1, ie+1
             if (itmask(i)) then
 c     Only calculate if H2I(i) is a substantial fraction
-              if (d(i,j,k)*dom.gt.1e10_RKIND) then
+              if (d(i,j,k)*dom > 1e10_RKIND) then
                 ciefudge = 1._RKIND
-                tau = ((d(i,j,k)/2e16_RKIND)*dom)**2.8_RKIND  ! 2e16 is in units of cm^-3
+                tau = ((d(i,j,k)/2.e16_RKIND)*dom)**2.8_RKIND  ! 2e16 is in units of cm^-3
                 tau = max(tau, 1.e-5_RKIND)
                 ciefudge = min((1._RKIND-exp(-tau))/tau,1._RKIND)
 c               Matt's attempt at a second exponentialier cutoff
 c     .                     (hdlte1/(1.0 + hdlte1/hdlow1)/(2.0*dom))
 c  new (correct) way: (april 4, 2007)
                hdlte1 = hdlte(i)/(HI(i,j,k)*dom)
-               hdlow1 = max(hdlow(i), REAL(tiny,RKIND))
+               hdlow1 = max(hdlow(i), tiny)
                edot(i) = edot(i) - HDI(i,j,k)*
      .              (hdlte(i)/(1._RKIND + hdlte1/hdlow1))/(3._RKIND*dom)
             end if
 
          logxe0 = log(xe_start)
          logxe1 = log(xe_end)
-         dlogxea = (logxe1 - logxe0) / REAL(n_xe-1,RKIND)
+         dlogxea = (logxe1 - logxe0) / REAL(n_xe-1)
 c
 c  Determine temperature for CMB temperature so we can subtract metal
 c  cooling at T=T_cmb
 c
          logtcmb = log(comp2)
-         tcmb_idx = max(1_IKIND, 
-     $        int((logtcmb-logtem0)/dlogtem,IKIND)+1_IKIND)
+         tcmb_idx = max(1, int((logtcmb-logtem0)/dlogtem,IKIND)+1)
          dlogTcmb = (logtcmb - (logtem0 + (tcmb_idx-1)*dlogtem)) 
      $        / dlogtem
          do i = is+1, ie+1
             if (itmask(i)) then
 !
-               if (de(i,j,k) .lt. 0._RKIND) then !#####
+               if (de(i,j,k) .lt. 0.0) then !#####
                   write(6,*) 'cool1d_multi: negative d(i,j,k)', 
      $                 de(i,j,k), d(i,j,k), metal(i,j,k)
                endif
      $                         xe_min)
                   log_xe = min(log_xe, xe_max)
                endif
-               xe_idx = min(n_xe-1, max(1_IKIND, 
-     $              int((log_xe - logxe0)/dlogxea,IKIND)+1_IKIND))
+               xe_idx = min(n_xe-1, max(1, 
+     $              int((log_xe - logxe0)/dlogxea,IKIND)+1))
                xe1 = logxe0 + (xe_idx - 1)*dlogxea
                dlogT = tdef(i)
                dlogxe = (log_xe - xe1) / dlogxea
                edot(i) = edot(i) - metalc(i) * d(i,j,k) * d(i,j,k) * 
      $              (xi/JHW_METALS_NORMZ)
 #ifdef UNUSED
-               if (xi.gt.1.e-4_RKIND .and. tgas(i).lt.1.e4_RKIND) then
+               if (xi > 1.e-4_RKIND .and. tgas(i) < 1.e4_RKIND) then
                   write(6,'(a,3i5,2e12.4)') 'a0', i,j,k, exp(log_xe),
      $                 exp(logtem(i))
                   write(6,'(a,4e12.4)') 'a1', d(i,j,k), xi, 
      $                        de(i,j,k), metal(i,j,k), d(i,j,k),  
      $                        HII(i,j,k), HeII(i,j,k), HeIII(i,j,k),
      $                        xe_max, xe_min
-                  edot(i) = 0._RKIND
+                  edot(i) = 0.0
                endif
 
             endif
          do j1=1, NID
             denwk(j1) = 10._RKIND**(DENMIN+(j1-0.5_RKIND)*DELD)
             do i=1, NIT
-               cbovcool(i,j1) = 0._RKIND
-               cbovheat(i,j1) = 0._RKIND
+               cbovcool(i,j1) = 0.0
+               cbovheat(i,j1) = 0.0
             enddo
          enddo
 
          else
 
             do i=1, NIB-1
-               radt(i) = 0._RKIND
+               radt(i) = 0.0
             enddo
 
          endif
 
-         radt(NIB) = 0._RKIND
+         radt(NIB) = 0.0
 
 !!! JB 2008-05-20 commented to avoid portability issues with Fortran MPI
 !!!         call open_mpi_error_file('Metal', 57, 'unknown')
 
 c           MKRJ: For low T, adopt the x=10^-2 cooling curve from
 c           Dalgarno & McCray 1972
-	    if (tgas(i) .ge. 1.e3) then
-               if (tgas(i) .ge. 6.31e3) then
-                  metal_cool = 3.13e-27*(tgas(i)**0.396_RKIND)
+	    if (tgas(i) .ge. 1.e3_RKIND) then
+               if (tgas(i) .ge. 6.31e3_RKIND) then
+                  metal_cool = 3.13e-27_RKIND*(tgas(i)**0.396_RKIND)
                else
-                  if (tgas(i) .ge. 3.16e3) then
-                     metal_cool = 2.64e-26*(tgas(i)**0.152_RKIND)
+                  if (tgas(i) .ge. 3.16e3_RKIND) then
+                     metal_cool = 2.64e-26_RKIND*(tgas(i)**0.152_RKIND)
                   else
-                     metal_cool = 5.28e-27*(tgas(i)**0.352_RKIND)
+                     metal_cool = 5.28e-27_RKIND*(tgas(i)**0.352_RKIND)
                   endif
                endif
             else
-               if (tgas(i) .ge. 39.8_RKIND) then
-                  if (tgas(i) .ge. 200._RKIND) then
-                     metal_cool = 3.06e-27*(tgas(i)**0.431_RKIND)
+               if (tgas(i) .ge. 3.98e1_RKIND) then
+                  if (tgas(i) .ge. 2.e2_RKIND) then
+                     metal_cool = 3.06e-27_RKIND*(tgas(i)**0.431_RKIND)
                   else
-                     metal_cool = 1.52e-28*(tgas(i)**0.997_RKIND)
+                     metal_cool = 1.52e-28_RKIND*(tgas(i)**0.997_RKIND)
                   endif
                else
-                  if (tgas(i) .ge. 25.1_RKIND) then
-                     metal_cool = 2.39e-29*(tgas(i)**1.5_RKIND)
+                  if (tgas(i) .ge. 2.51e1_RKIND) then
+                     metal_cool = 2.39e-29_RKIND*(tgas(i)**1.5_RKIND)
                   else
-                     metal_cool = 1.095e-32*(tgas(i)**3.885_RKIND)
+                     metal_cool = 1.095e-32_RKIND*(tgas(i)**3.885_RKIND)
                   endif
                endif
             endif

File src/enzo/cool1d_sep.F

 c
 c  Parameters
 c
-      double precision mh
+      real*8 mh
       R_PREC    ZSOLAR
       parameter (mh = mass_h, ZSOLAR = 0.02041_RKIND)
 c
       R_PREC dt2, ttmin, comp1, comp2, scoef, acoef, HIdot,
      &     hdlte1, hdlow1, gamma2, x, nH2, nother, fudge, fH2,
      &     gphdl1, factor, ciefudge, tau, dlogTcmb, logtcmb, metalc_cmb
-      double precision coolunit, dbase1, tbase1, xbase1, rtunits
+      real*8 coolunit, dbase1, tbase1, xbase1, rtunits
 c
 c  Slice locals
 c 
       INTG_PREC indixe(in)
       R_PREC t1(in), t2(in), logtem(in), tdef(in), p2d(in),
      &     tgas(in), tgasold(in)
-      double precision edot(in,nedot)
+      real*8 edot(in,nedot)
 c
 c  Cooling/heating slice locals
 c
-      double precision ceHI(in), ceHeI(in), ceHeII(in),
+      real*8 ceHI(in), ceHeI(in), ceHeII(in),
      &     ciHI(in), ciHeI(in), ciHeIS(in), ciHeII(in),
      &     reHII(in), reHeII1(in), reHeII2(in), reHeIII(in),
      &     brem(in)
 c
 c        Compute index into the table and precompute parts of linear interp
 c
-         indixe(i) = min(nratec-1, max(1_IKIND,
-     &        int((logtem(i)-logtem0)/dlogtem,IKIND)+1_IKIND))
+         indixe(i) = min(nratec-1, max(1,
+     &        int((logtem(i)-logtem0)/dlogtem,IKIND)+1))
          t1(i) = (logtem0 + (indixe(i) - 1)*dlogtem)
          t2(i) = (logtem0 + (indixe(i)    )*dlogtem)
          tdef(i) = 1._RKIND/(t2(i) - t1(i))
 c  cooling at T=T_cmb
 c
          logtcmb = log(comp2)
-         tcmb_idx = max(1_IKIND, 
-     $        int((logtcmb-logtem0)/dlogtem,IKIND)+1_IKIND)
+         tcmb_idx = max(1, int((logtcmb-logtem0)/dlogtem,IKIND)+1)
          dlogTcmb = (logtcmb - (logtem0 + (tcmb_idx-1)*dlogtem)) 
      $	/ dlogtem
          do i = is+1, ie+1
      $                         xe_min)
                   log_xe = min(log_xe, xe_max)
                endif
-               xe_idx = min(n_xe-1, max(1_IKIND, 
-     $              int((log_xe - logxe0)/dlogxea,IKIND)+1_IKIND))
+               xe_idx = min(n_xe-1, max(1, 
+     $              int((log_xe - logxe0)/dlogxea,IKIND)+1))
                xe1 = logxe0 + (xe_idx - 1)*dlogxea
                dlogT = (logtem(i) - t1(i)) * tdef(i)
                dlogxe = (log_xe - xe1) / dlogxea

File src/enzo/cool_multi_lum.F

 c
       INTG_PREC i, j, k, n
       R_PREC comp1, comp2, energy
-      double precision dom, coolunit, dbase1, tbase1, xbase1, mh, ulum
+      real*8 dom, coolunit, dbase1, tbase1, xbase1, mh, ulum
       parameter (mh = mass_h)
 c
 c  Slice locals
       INTG_PREC indixe(ijk)
       R_PREC t1(ijk), t2(ijk), logtem(ijk), tdef(ijk), p2d(ijk),
      &     dtit(ijk), ttot(ijk), tgas(ijk), tgasold(ijk)
-      double precision edot(in,nlum)
+      real*8 edot(in,nlum)
 c
 c  Cooling/heating slice locals
 c
-      double precision ceHI(ijk), ceHeI(ijk), ceHeII(ijk),
+      real*8 ceHI(ijk), ceHeI(ijk), ceHeII(ijk),
      &     ciHI(ijk), ciHeI(ijk), ciHeIS(ijk), ciHeII(ijk),
      &     reHII(ijk), reHeII1(ijk), reHeII2(ijk), reHeIII(ijk),
      &     brem(ijk)

File src/enzo/cool_multi_time.F

       R_PREC t1(ijk), t2(ijk), logtem(ijk), tdef(ijk), p2d(ijk),
      &     dtit(ijk), ttot(ijk), tgas(ijk), tgasold(ijk),
      &     tdust(ijk), metallicity(ijk), rhoH(ijk)
-      double precision edot(ijk)
+      real*8 edot(ijk)
 
 !  Cooling/heating slice locals
 
-      double precision ceHI(ijk), ceHeI(ijk), ceHeII(ijk),
+      real*8 ceHI(ijk), ceHeI(ijk), ceHeII(ijk),
      &     ciHI(ijk), ciHeI(ijk), ciHeIS(ijk), ciHeII(ijk),
      &     reHII(ijk), reHeII1(ijk), reHeII2(ijk), reHeIII(ijk),
      &     brem(ijk), cieco(ijk)

File src/enzo/fourn.F

      &        ifp1, ifp2, ip1, ip2, ip3, k2, n, nprev, nrem, ntot
       INTG_PREC ibad
       R_PREC    tempi, tempr
-c     double precision twopi, theta, wr, wi, wpr, wpi, wtemp
+c     real*8 twopi, theta, wr, wi, wpr, wpi, wtemp
       R_PREC    twopi, theta, wr, wi, wpr, wpi, wtemp
       parameter (twopi=6.2831853071796_RKIND)
 c

File src/enzo/gFLDProblem_AnalyticalEqns.F90

      lTempE = log(TempEnd)
      dlTemp = (lTempE - lTempS)/(1.d0*NTempBins - 1.d0)
      lTemp = min(max(log(T), lTempS), lTempE)
-     Tidx = min(NTempBins-1, max(1_IKIND, int((lTemp-lTempS)/dlTemp,IKIND)+1))
+     Tidx = min(NTempBins-1, max(1, int((lTemp-lTempS)/dlTemp,IKIND)+1))
      Tidxp = Tidx+1
      Tl = lTempS + (Tidx-1)*dlTemp
      Tr = lTempS +  Tidx*dlTemp
      lTempE = log(TempEnd)
      dlTemp = (lTempE - lTempS)/(1.d0*NTempBins - 1.d0)
      lTemp = min(max(log(T), lTempS), lTempE)
-     Tidx = min(NTempBins-1, max(1_IKIND, int((lTemp-lTempS)/dlTemp,IKIND)+1))
+     Tidx = min(NTempBins-1, max(1, int((lTemp-lTempS)/dlTemp,IKIND)+1))
      Tidxp = Tidx+1
      Tl = lTempS + (Tidx-1)*dlTemp
      Tr = lTempS +  Tidx*dlTemp
      T = max(T,min_temp)
      lamT = 3.15614e5_RKIND/T
      lTemp = min(max(log(T), lTempS), lTempE)
-     Tidx = min(NTempBins-1, max(1_IKIND, int((lTemp-lTempS)/dlTemp,IKIND)+1))
+     Tidx = min(NTempBins-1, max(1, int((lTemp-lTempS)/dlTemp,IKIND)+1))
      Tidxp = Tidx+1
      Tl = lTempS + (Tidx-1)*dlTemp
      Tr = lTempS +  Tidx*dlTemp
      T = max(T,min_temp)
      lamT = 3.15614d5/T
      lTemp = min(max(log(T), lTempS), lTempE)
-     Tidx = min(NTempBins-1, max(1_IKIND, int((lTemp-lTempS)/dlTemp,IKIND)+1))
+     Tidx = min(NTempBins-1, max(1, int((lTemp-lTempS)/dlTemp,IKIND)+1))
      Tidxp = Tidx+1
      Tl = lTempS + (Tidx-1)*dlTemp
      Tr = lTempS +  Tidx*dlTemp

File src/enzo/gFLDProblem_LocalJac.F90

 
               ! look up rates, derivatives
               lTemp = min(max(log(T), lTempS), lTempE)
-              Tidx = min(NTempBins-1, max(1_IKIND, int((lTemp-lTempS)/dlTemp,IKIND)+1_IKIND))
+              Tidx = min(NTempBins-1, max(1, int((lTemp-lTempS)/dlTemp,IKIND)+1))
               Tidxp = Tidx+1
               Tl = lTempS + (Tidx-1)*dlTemp
               Tr = lTempS +  Tidx*dlTemp
 
               ! look up rates
               lTemp = min(max(log(T), lTempS), lTempE)
-              Tidx = min(NTempBins-1, max(1_IKIND, int((lTemp-lTempS)/dlTemp,IKIND)+1_IKIND))
+              Tidx = min(NTempBins-1, max(1, int((lTemp-lTempS)/dlTemp,IKIND)+1))
               Tidxp = Tidx+1
               Tl = lTempS + (Tidx-1)*dlTemp
               Tr = lTempS +  Tidx*dlTemp
 
               ! look up rates, derivatives
               lTemp = min(max(log(T), lTempS), lTempE)
-              Tidx = min(NTempBins-1, max(1_IKIND, int((lTemp-lTempS)/dlTemp,IKIND)+1_IKIND))
+              Tidx = min(NTempBins-1, max(1, int((lTemp-lTempS)/dlTemp,IKIND)+1))
               Tidxp = Tidx+1
               Tl = lTempS + (Tidx-1)*dlTemp
               Tr = lTempS +  Tidx*dlTemp
 
               ! look up rates
               lTemp = min(max(log(T), lTempS), lTempE)
-              Tidx = min(NTempBins-1, max(1_IKIND, int((lTemp-lTempS)/dlTemp,IKIND)+1_IKIND))
+              Tidx = min(NTempBins-1, max(1, int((lTemp-lTempS)/dlTemp,IKIND)+1))
               Tidxp = Tidx+1
               Tl = lTempS + (Tidx-1)*dlTemp
               Tr = lTempS +  Tidx*dlTemp
 
            ! look up rates
            lTemp = min(max(log(T), lTempS), lTempE)
-           Tidx = min(NTempBins-1, max(1_IKIND, int((lTemp-lTempS)/dlTemp,IKIND)+1_IKIND))
+           Tidx = min(NTempBins-1, max(1, int((lTemp-lTempS)/dlTemp,IKIND)+1))
            Tidxp = Tidx+1
            Tl = lTempS + (Tidx-1)*dlTemp
            Tr = lTempS +  Tidx*dlTemp

File src/enzo/gFLDProblem_LocalRHS.F90

 
               ! look up rates
               lTemp = min(max(log(T), lTempS), lTempE)
-              Tidx = min(NTempBins-1, max(1_IKIND, int((lTemp-lTempS)/dlTemp,IKIND)+1_IKIND))
+              Tidx = min(NTempBins-1, max(1, int((lTemp-lTempS)/dlTemp,IKIND)+1))
               Tidxp = Tidx+1
               Tl = lTempS + (Tidx-1)*dlTemp
               Tr = lTempS +  Tidx*dlTemp
                  
                  ! look up rates
                  lTemp = min(max(log(T), lTempS), lTempE)
-                 Tidx = min(NTempBins-1, max(1_IKIND, int((lTemp-lTempS)/dlTemp,IKIND)+1_IKIND))
+                 Tidx = min(NTempBins-1, max(1, int((lTemp-lTempS)/dlTemp,IKIND)+1))
                  Tidxp = Tidx+1
                  Tl = lTempS + (Tidx-1)*dlTemp
                  Tr = lTempS +  Tidx*dlTemp
                  
                  ! look up rates
                  lTemp = min(max(log(T), lTempS), lTempE)
-                 Tidx = min(NTempBins-1, max(1_IKIND, int((lTemp-lTempS)/dlTemp,IKIND)+1_IKIND))
+                 Tidx = min(NTempBins-1, max(1, int((lTemp-lTempS)/dlTemp,IKIND)+1))
                  Tidxp = Tidx+1
                  Tl = lTempS + (Tidx-1)*dlTemp
                  Tr = lTempS +  Tidx*dlTemp
                  
                  ! look up rates
                  lTemp = min(max(log(T), lTempS), lTempE)
-                 Tidx = min(NTempBins-1, max(1_IKIND, int((lTemp-lTempS)/dlTemp,IKIND)+1_IKIND))
+                 Tidx = min(NTempBins-1, max(1, int((lTemp-lTempS)/dlTemp,IKIND)+1))
                  Tidxp = Tidx+1
                  Tl = lTempS + (Tidx-1)*dlTemp
                  Tr = lTempS +  Tidx*dlTemp
                  
                  ! look up rates
                  lTemp = min(max(log(T), lTempS), lTempE)
-                 Tidx = min(NTempBins-1, max(1_IKIND, int((lTemp-lTempS)/dlTemp,IKIND)+1_IKIND))
+                 Tidx = min(NTempBins-1, max(1, int((lTemp-lTempS)/dlTemp,IKIND)+1))
                  Tidxp = Tidx+1
                  Tl = lTempS + (Tidx-1)*dlTemp
                  Tr = lTempS +  Tidx*dlTemp

File src/enzo/gFLDSplit_AnalyticChemistry.F90

   lTempE = log(TempEnd)
   dlTemp = (lTempE - lTempS)/(1.d0*NTempBins - 1.d0)
   lTemp = min(max(log(T), lTempS), lTempE)
-  Tidx = min(NTempBins-1, max(1_IKIND, int((lTemp-lTempS)/dlTemp,IKIND)+1_IKIND))
+  Tidx = min(NTempBins-1, max(1, int((lTemp-lTempS)/dlTemp,IKIND)+1))
   Tidxp = Tidx+1
   Tl = lTempS + (Tidx-1)*dlTemp
   Tr = lTempS +  Tidx*dlTemp
      lTempE = log(TempEnd)
      dlTemp = (lTempE - lTempS)/(1.d0*NTempBins - 1.d0)
      lTemp = min(max(log(T), lTempS), lTempE)
-     Tidx = min(NTempBins-1, max(1_IKIND, int((lTemp-lTempS)/dlTemp,IKIND)+1_IKIND))
+     Tidx = min(NTempBins-1, max(1, int((lTemp-lTempS)/dlTemp,IKIND)+1))
      Tidxp = Tidx+1
      Tl = lTempS + (Tidx-1)*dlTemp
      Tr = lTempS +  Tidx*dlTemp

File src/enzo/interpolate.F

 c  Error check for self-consistent interpolation indexes
 c
       do i = 1, ndim
-         if (int(pstart(i)/refine(i))*refine(i) .ne. pstart(i)) then
+         if (int(pstart(i)/refine(i),IKIND)*refine(i) /= pstart(i)) then
             write(6,*) "INTERPOLATE: ERROR pstart(i) =", pstart(i),
      &                 " refine(i) =", refine(i), " i =", i
             ERROR_MESSAGE

File src/enzo/mcooling.F

 #include "phys_const.def"
 #include "error.def"
 
+#ifdef CEN_METALS
 C     NAME MTLCLHT
       SUBROUTINE MTLCLHT(DEN,RADT,IPHOTT,CBOVCOOL,CBOVHEAT)
-#ifdef CEN_METALS
 C
 C     THIS PROGRAM CALCULATES A COOLING/HEATING TABLE FOR METALS FROM
 C     C TO NI EXCLUDING H AND HE, GIVEN AN IONIZING RADIATION FIELD
 C     R. MOORE OCTOBER 1976
 #include "fortran_types.def"
       IMPLICIT NONE
-      DOUBLE PRECISION XD, XS1, XS2, X3, A10, A11, A12, A13, A14, A15, 
+      REAL*8 XD, XS1, XS2, X3, A10, A11, A12, A13, A14, A15, 
      &                 A20, A21, A22, A23, A24, A25, A30, A31, A32, A33, 
      &                 A34, A35, B10, B11, B12, B13, B14, B20, B21, B22, 
      &                 B23, B24, B30, B31, B32, B33, B34, SEATON
      .                      +(32.d0-448.d0/X3)/X3)/X3)
    50 X3=X**0.33333333333333d0*Q**0.66666666666667d0
       SEATON=EXINT1(X,3)+(XS1+XS2/X3)/X3
-#endif /* CEN_METALS */
       RETURN
       END
+#endif /* CEN_METALS */

File src/enzo/mg_restrict.F

          enddo
          do i=1, ddim1
             i1 = min(max(int(REAL(i-1,RKIND)*fact1 + 0.5_RKIND,IKIND)+1, 
-     &           1_IKIND), sdim1)
+     &           1), sdim1)
             dest(i,    1,1) = source(i1,    1,1)
             dest(i,ddim2,1) = source(i1,sdim2,1)
          enddo
             enddo
             do i=1, ddim1
                i1 = min(max(int(REAL(i-1,RKIND)*fact1 + 0.5_RKIND,IKIND) 
-     &              + 1, 1_IKIND), sdim1)
+     &              + 1, 1), sdim1)
                dest(i,    1,k) = source(i1,    1,k1)
                dest(i,ddim2,k) = source(i1,sdim2,k1)
             enddo
          enddo
          do j=1, ddim2
             j1 = min(max(int(REAL(j-1,RKIND)*fact2 + 0.5_RKIND,IKIND)+1, 
-     &           1_IKIND), sdim2)
+     &           1), sdim2)
             do i=1, ddim1
                i1 = min(max(int(REAL(i-1,RKIND)*fact1 + 0.5_RKIND,IKIND) 
-     &              + 1, 1_IKIND), sdim1)
+     &              + 1, 1), sdim1)
                dest(i,j,    1) = source(i1,j1,    1)
                dest(i,j,ddim3) = source(i1,j1,sdim3)
             enddo

File src/enzo/multi_cool.F

       parameter (tolerance = 1.e-10_RKIND)
 #endif
 
-      double precision mh
+      real*8 mh
       parameter (mh = mass_h)
 
 !  Locals
       R_PREC t1(ijk), t2(ijk), logtem(ijk), tdef(ijk), p2d(ijk),
      &     dtit(ijk), ttot(ijk), tgas(ijk), tgasold(ijk),
      &     tdust(ijk), metallicity(ijk), rhoH(ijk)
-      double precision edot(ijk)
+      real*8 edot(ijk)
 
 !  Cooling/heating row locals
 
-      double precision ceHI(ijk), ceHeI(ijk), ceHeII(ijk),
+      real*8 ceHI(ijk), ceHeI(ijk), ceHeII(ijk),
      &     ciHI(ijk), ciHeI(ijk), ciHeIS(ijk), ciHeII(ijk),
      &     reHII(ijk), reHeII1(ijk), reHeII2(ijk), reHeIII(ijk),
      &     brem(ijk), cieco(ijk)

File src/enzo/pop3_color_maker.F

       R_PREC   div, tdyn, dtot, wgt
       R_PREC   pi, G, sndspdC, m1
       R_PREC   isosndsp2, bmass, jeanmass, lum
-      double precision msolar, mh
+      real*8 msolar, mh
       parameter (pi=pi_val, G=GravConst, 
      &           sndspdC=1.3095e8_RKIND,
      &           msolar=SolarMass, mh=mass_h)

File src/enzo/pop3_maker.F

       R_PREC   div, tdyn, dtot, wgt
       R_PREC   pi, G, sndspdC, m1
       R_PREC   isosndsp2, bmass, jeanmass, lum
-      double precision msolar, mh
+      real*8 msolar, mh
       parameter (pi=pi_val, G=GravConst, 
      &           sndspdC=1.3095e8_RKIND,
      &           msolar=SolarMass, mh=mass_h)
       INTG_PREC i, j, k, n
       R_PREC energy, lum, sn_nrg, m_eject, yield, dratio, m1
       R_PREC lower_pisn, upper_pisn, lsun, he_core, starmass, lifetime
-      double precision msolar
+      real*8 msolar
       parameter (lower_pisn = 140._RKIND, upper_pisn = 260._RKIND, 
      $           lsun = 2.e33_RKIND)
       parameter (msolar = SolarMass)

File src/enzo/rotate2d.F

       INTG_PREC bs
 
       bs = 64
-      ib = n1/max(n1/bs,1_IKIND)
-      jb = n2/max(n2/bs,1_IKIND)
+      ib = n1/max(n1/bs,1)
+      jb = n2/max(n2/bs,1)
 
 !     do i=1,n1
 !     do j=1,n2

File src/enzo/rotate3d.F

 
       bs = 64
 
-      ib=n1/max(n1/bs,1_IKIND)
-      jb=n2/max(n2/bs,1_IKIND)
-      kb=n3/max(n3/bs,1_IKIND)
+      ib=n1/max(n1/bs,1)
+      jb=n2/max(n2/bs,1)
+      kb=n3/max(n3/bs,1)
 
 !     do i=1,n1
 !     do k=1,n3

File src/enzo/smooth.F

                enddo
                do n=-nsmooth,+nsmooth
                   do i=1, sdim1
-                     i1 = min(max(i+n, 2_IKIND), sdim1)
+                     i1 = min(max(i+n, 2), sdim1)
                      temp(i) = temp(i) + source1(i1,j,k)
                   enddo
                enddo
                enddo
                do n=-nsmooth,+nsmooth
                   do j=1, sdim2
-                     j1 = min(max(j+n, 2_IKIND), sdim2)
+                     j1 = min(max(j+n, 2), sdim2)
                      temp(j) = temp(j) + source2(i,j1,k)
                   enddo
                enddo
                enddo
                do n=-nsmooth,+nsmooth
                   do k=1, sdim3
-                     k1 = min(max(k+n, 2_IKIND), sdim3)
+                     k1 = min(max(k+n, 2), sdim3)
                      temp(k) = temp(k) + source3(i,j,k1)
                   enddo
                enddo

File src/enzo/solve_rate.F

       parameter (tolerance = 1.e-10_RKIND)
 #endif
 
-      double precision everg, mh, e24, e26
+      real*8 everg, mh, e24, e26
       parameter (everg = ev2erg, mh = mass_h,
      &           e24 = 13.6_RKIND, e26 = 24.6_RKIND)            ! DPC
 c
      &     totalH(ijk), totalHe(ijk), totalD, 
      &     correctH, correctHe, correctD, metalfree
 
-      double precision coolunit, dbase1, tbase1, xbase1
+      real*8 coolunit, dbase1, tbase1, xbase1
 c
 c  Slice temporaries (passed in)
 c 
 c
 c           Find index into tble and precompute interpolation values
 c
-            indixe(i) = min(nratec-1, max(1_IKIND,
-     &           int((logtem(i)-logtem0)/dlogtem,IKIND)+1_IKIND))
+            indixe(i) = min(nratec-1, max(1,
+     &           int((logtem(i)-logtem0)/dlogtem,IKIND)+1))
             t1(i) = (logtem0 + (indixe(i) - 1)*dlogtem)
             t2(i) = (logtem0 + (indixe(i)    )*dlogtem)
             tdef(i) = t2(i) - t1(i)

File src/enzo/solve_rate_cool.F

       parameter (tolerance = 1.e-10_RKIND)
 #endif
 
-      double precision mh
+      real*8 mh
       parameter (mh = mass_h)
 
 !  Locals
       INTG_PREC i, j, k, iter
       INTG_PREC clGridDim1, clGridDim2, clGridDim3, clGridDim4, 
      &     clGridDim5
-      R_PREC ttmin, dom, energy, comp1, comp2, dlogtem
-      double precision coolunit, dbase1, tbase1, xbase1, chunit, uvel
-      double precision heq1, heq2, eqk221, eqk222, eqk131, eqk132,
-     &                 eqt1, eqt2, eqtdef, dheq, heq
+      R_PREC ttmin, dom, energy, comp1, comp2
+      real*8 coolunit, dbase1, tbase1, xbase1, chunit, uvel
+      real*8 heq1, heq2, eqk221, eqk222, eqk131, eqk132,
+     &                 eqt1, eqt2, eqtdef, dheq, heq, dlogtem
 
 !  row temporaries
 
 
 !  Cooling/heating row locals
 
-      double precision ceHI(ijk), ceHeI(ijk), ceHeII(ijk),
+      real*8 ceHI(ijk), ceHeI(ijk), ceHeII(ijk),
      &     ciHI(ijk), ciHeI(ijk), ciHeIS(ijk), ciHeII(ijk),
-     &     reHII(ijk), reHeII1(ijk), reHeII2(ijk), reHeIII(ijk),
+     &     reHII(ijk), reHeII1(ijk), reHeII2(ijk), reHeIII(ijk), 
      &     brem(ijk), edot(ijk)
-      R_PREC hy_RKIND1k(ijk), h2k01(ijk), vibh(ijk), roth(ijk), 
-     &     rotl(ijk), gpldl(ijk), gphdl(ijk), hdlte(ijk), hdlow(ijk), 
-     &     cieco(ijk)
+      R_PREC hyd01k(ijk), h2k01(ijk), vibh(ijk), roth(ijk), rotl(ijk), 
+     &     gpldl(ijk), gphdl(ijk), hdlte(ijk), hdlow(ijk), cieco(ijk)
 
 !  Iteration mask
 
      &                piHI, piHeI, piHeII, comp1, comp2,
      &                HM, H2I, H2II, DI, DII, HDI, metal,
      &                hyd01ka, h2k01a, vibha, rotha, rotla,
-     &                hy_RKIND1k, h2k01, vibh, roth, rotl,
+     &                hyd01k, h2k01, vibh, roth, rotl,
      &                gpldla, gphdla, gpldl, gphdl,
      &                hdltea, hdlowa, hdlte, hdlow,
      &                gaHIa, gaH2a, gaHea, gaHpa, gaela,
      &        idust, ndratec
       R_PREC temstart, temend, tgas1d(in), dom,
      &        dtemstart, dtemend
-      double precision coolunit, tbase1
+      real*8 coolunit, tbase1
       logical itmask(ijk)
 
 !     Chemistry rates as a function of temperature
 
 !     Parameters
 
-      double precision everg, e24, e26
+      real*8 everg, e24, e26
       parameter(everg = ev2erg, e24 = 13.6_RKIND, 
      &          e26 = 24.6_RKIND)
 
 
 !        Find index into tble and precompute interpolation values
 
-         indixe(i) = min(nratec-1, max(1_IKIND,
+         indixe(i) = min(nratec-1, max(1,
      &        int((logtem(i)-logtem0)/dlogtem,IKIND)+1))
          t1(i) = (logtem0 + (indixe(i) - 1)*dlogtem)
          t2(i) = (logtem0 + (indixe(i)    )*dlogtem)
 
 !                 Find index into table and precompute interpolation values
 
-                  d_indixe(i) = min(ndratec-1,
-     &                 max(1,int((d_logtem(i)-d_logtem0)/d_dlogtem)+1))
+                  d_indixe(i) = min(ndratec-1, max(1,
+     &                 int((d_logtem(i)-d_logtem0)/d_dlogtem,IKIND)+1))
                   d_t1(i) = (d_logtem0 + (d_indixe(i) - 1)*d_dlogtem)
                   d_t2(i) = (d_logtem0 + (d_indixe(i)    )*d_dlogtem)
                   d_tdef(i) = (d_logtem(i) - d_t1(i)) / 
       INTG_PREC ispecies, idust, is, ie, j, k, in, jn, kn, 
      &     iradtrans, irt_honly
       R_PREC dedot(in), HIdot(in), dom
-      double precision edot(in)
+      real*8 edot(in)
       logical itmask(in)
 
 !     Density fields
       R_PREC    kphHI(in,jn,kn), kphHeI(in,jn,kn), kphHeII(in,jn,kn),
      &        kdissH2I(in,jn,kn)
 
-      double precision chunit
+      real*8 chunit
 
 !     Rate values
 
 !     locals
 
       INTG_PREC i
-      double precision h2heatfac(in), H2delta(in), H2dmag, atten, tau
+      real*8 h2heatfac(in), H2delta(in), H2dmag, atten, tau
 
       if (ispecies .eq. 1) then
 

File src/enzo/star_feedback_pn_snia.F

       R_PREC MdotSNII, MdotSNIa, MdotPN
       R_PREC Edota, Edotb, Mdota, Mdotb, dvol
       R_PREC tau, rmf, dtt, madd, madd_cell, delta_e, tmid, tnow
-      double precision msolar, esnIa, eSNII
+      real*8 msolar, esnIa, eSNII
       parameter (tau=1.0e9_RKIND,  msolar=1.989D33) ! tau in years
       parameter (eSNIa=5.D-6, eSNII=5.D-6) 
 c-----------------------------------------------------------------------
 c
 c        Compute index of the cell that the star particle resides in.
 c 
-         i = int((xp(n) - xstart)/dx) + 1
-         j = int((yp(n) - ystart)/dx) + 1
-         k = int((zp(n) - zstart)/dx) + 1
+         i = int((xp(n) - xstart)/dx,IKIND) + 1
+         j = int((yp(n) - ystart)/dx,IKIND) + 1
+         k = int((zp(n) - zstart)/dx,IKIND) + 1
 c
 c         check bounds - if star particle is outside of this grid
 c         then exit and give a warning.

File src/enzo/star_maker1.F

       R_PREC   div, tdyn, dtot
       R_PREC   pi, G, sndspdC
       R_PREC   isosndsp2, starmass, starfraction, bmass, jeanmass
-      double precision msolar
+      real*8 msolar
       parameter (pi=pi_val, G=GravConst, 
      &           sndspdC=1.3095e8_RKIND,
      &           msolar=SolarMass)

File src/enzo/star_maker2.F

       R_PREC   div, tdyn, dtot
       R_PREC   pi, G, sndspdC
       R_PREC   isosndsp2, starmass, starfraction, bmass, jeanmass
-      double precision msolar
+      real*8 msolar
       parameter (pi=pi_val, G=GravConst, 
      &           sndspdC=1.3095e8_RKIND,
      &           msolar=SolarMass)

File src/enzo/star_maker3.F

       R_PREC   div, tdyn, dtot
       R_PREC   pi, G, sndspdC
       R_PREC   isosndsp2, starmass, starfraction, bmass, jeanmass
-      double precision msolar
+      real*8 msolar
       parameter (pi=pi_val, G=GravConst, 
      &           sndspdC=1.3095e8_RKIND,
      &           msolar=SolarMass)

File src/enzo/star_maker4.F

       R_PREC   pi, G, sndspdC,mproton
       R_PREC   isosndsp2, starmass, starfraction, bmass, jeanmass
       R_PREC   densthresh, timeconstant, gasfrac
-      double precision msolar
+      real*8 msolar
       parameter (pi=pi_val, G=GravConst, 
      &           sndspdC=1.3095e8_RKIND,
      &           msolar=SolarMass,

File src/enzo/star_maker5.F

       R_PREC   usn, tstar, pstar, sqrtepssn, timeconstant
       R_PREC   mproton, RandomSeed, beta,x ,y, kb
       R_PREC   starfraction, bmass, densthresh 
-      double precision msolar
+      real*8 msolar
       parameter (msolar=SolarMass, mproton=mass_h, 
      &           beta=0.1_RKIND, sqrtepssn=2.e24, 
      &           RandomSeed=-123456, kb=kboltz)

File src/enzo/star_maker7.F

       R_PREC   div, tdyn, dtot
       R_PREC   pi, G, sndspdC
       R_PREC   isosndsp2, starmass, starfraction, bmass, jeanmass
-      double precision msolar
+      real*8 msolar
       parameter (pi=pi_val, G=GravConst, sndspdC=1.3095d8,
      &           msolar=SolarMass)
 c

File src/enzo/star_maker_h2reg.F

 
       R_PREC U_MW, D_MW, Dstar, alpha, gg, Lambda, x
 
-      double precision msolar
+      real*8 msolar
       parameter (pi=pi_val, G=6.67e-8_RKIND,
      $     msolar=1.989e33_RKIND,mproton=1.67e-24_RKIND)