fix precision loss in complex number comversion in turbulence module

Author: tinecolman

tests/turb/driving/driving-ref.dat

@@ -1,16 +1,16 @@
 boxlen : 4.0000000000000000e+00
-density : -1.5970192481745945e+05
+density : -1.5970192503110823e+05
 dx : 1.6384000000000000e+04
 level : 1.5728640000000000e+06
 ncells : 2.6214400000000000e+05
-pressure : -5.3599210383977520e+05
+pressure : -5.3599210405342397e+05
 time : 2.5147260806143501e-01
 unit_d : 1.5049295743500001e-20
 unit_l : 3.0857000000000000e+18
 unit_t : 3.1556926000000000e+13
-velocity_x : 2.2878501962849146e+05
-velocity_y : 2.2205591221498317e+05
-velocity_z : 2.1700854779837615e+05
+velocity_x : 2.2878502213449072e+05
+velocity_y : 2.2205591253462457e+05
+velocity_z : 2.1700854906616427e+05
 x : 5.2428800000000000e+05
 y : 5.2428800000000000e+05
 z : 5.2428800000000000e+05

turb/init_turb.f90

@@ -144,7 +144,7 @@ subroutine init_turb
       afield_next = afield_next * turb_norm * turb_rms
    else
       ! Not a restart - set up initial field and perform FFT
-      turb_next = cmplx(0.0_dp, 0.0_dp)
+      turb_next = cmplx(0, 0, kind=cdp)
       call add_turbulence(turb_next, turb_dt)
 
       ! Fourier transform

turb/turb_force_utils.f90

@@ -153,7 +153,7 @@ subroutine gaussian_cmplx(G)
    call kiss64_double(ndim, kiss64_state, Rnd(1:ndim))
    arg = (Rnd(1:ndim) * 2.0 * PI) - 1.0
 
-   G = cmplx(mag * cos(arg), mag * sin(arg))
+   G = cmplx(mag * cos(arg), mag * sin(arg), kind=cdp)
 
    ! Method 2: Gaussian real and imaginary components
 !    do d=1,ndim
@@ -187,19 +187,21 @@ subroutine add_turbulence(turb_field, dt)
                                            ! Wiener process
 
    integer              :: i, j, k         ! Loop variables
-   integer              :: ii,jj,kk        ! Hermitian pair
    integer              :: halfway         ! Half of grid size
 
 #if defined(HERMITIAN_FIELD)
+   integer              :: ii,jj,kk        ! Hermitian pair
    logical              :: hermitian_pair  ! Do we need to take conjugate?
    logical              :: own_conjg       ! Are we are own conjugate pair?
 #endif
 
-   real(kind=dp)        :: unitk(1:NDIM)       ! Unit vector parallel to k
    complex(kind=cdp)    :: complex_vec(1:NDIM) ! Complex Fourier vector
+#if NDIM>1
+   real(kind=dp)        :: unitk(1:NDIM)       ! Unit vector parallel to k
    complex(kind=cdp)    :: unitk_cmplx(1:NDIM) ! Unit vector parallel to k
    complex(kind=cdp)    :: comp_cmplx(1:NDIM)  ! Compressive component
    complex(kind=cdp)    :: sol_cmplx(1:NDIM)   ! Solenoidal components
+#endif
 
    halfway = TURB_GS / 2 ! integer division
 
@@ -262,7 +264,7 @@ subroutine add_turbulence(turb_field, dt)
 #if NDIM>1
             ! Helmholtz decomposition!
             call find_unitk(i, j, k, halfway, unitk)
-            unitk_cmplx = cmplx(unitk)
+            unitk_cmplx = cmplx(unitk, kind=cdp)
             comp_cmplx = unitk_cmplx *  dot_product(unitk_cmplx,complex_vec)
             sol_cmplx = complex_vec - comp_cmplx
             complex_vec = comp_cmplx*comp_frac + sol_cmplx*sol_frac
@@ -490,7 +492,7 @@ subroutine power_rms_norm(power_in, P)
    integer                    :: d               ! Dimension counter
 
    do d=1,NDIM
-      complex_field(d,:,:,:) = cmplx(power_in)
+      complex_field(d,:,:,:) = cmplx(power_in, kind=cdp)
    end do
 
 #if NDIM==1