Can't replace Fortran real variables by double precision variables or more precision

Question

I am using a known code (CAMB) which generates values like this :

k(h/Mpc) Pk/s8^2(Mpc/h)^3
5.2781500000e-06 1.9477400000e+01
5.5479700000e-06 2.0432300000e+01
5.8315700000e-06 2.1434000000e+01
6.1296700000e-06 2.2484700000e+01
6.4430100000e-06 2.3587000000e+01
6.7723700000e-06 2.4743400000e+01
7.1185600000e-06 2.5956400000e+01
7.4824500000e-06 2.7228900000e+01
7.8649500000e-06 2.8563800000e+01
8.2669900000e-06 2.9964100000e+01

I would like to get more precision on the generated values, like this :

k(h/Mpc) Pk/s8^2(Mpc/h)^3
5.3594794736e-06 1.8529569626e+01
5.6332442000e-06 1.9437295914e+01
5.9209928622e-06 2.0389484405e+01
6.2234403231e-06 2.1388326645e+01
6.5413364609e-06 2.2436098099e+01
6.8754711720e-06 2.3535198212e+01
7.2266739153e-06 2.4688137054e+01
7.5958159869e-06 2.5897554398e+01
7.9838137026e-06 2.7166225433e+01
8.3916311269e-06 2.8497039795e+01
8.8202796178e-06 2.9893053055e+01
9.2708232842e-06 3.1357446670e+01
9.7443817140e-06 3.2893573761e+01

Here the section of code that produces the data :

I tried to do the following modifications in the declarations of variables at the beginning of code above :

1)First try :

!Export files of total  matter power spectra in h^{-1} Mpc units, against k/h.
Type(MatterTransferData), intent(in) :: MTrans
Type(CAMBdata) :: State
character(LEN=Ini_max_string_len), intent(IN) :: FileNames(*)
character(LEN=name_tag_len) :: columns(3)
integer itf, i, unit
integer points
! Added : way of declaring double precision
integer, parameter :: wp = selected_real_kind(15,307)
real(wp), dimension(:,:), allocatable :: outpower

but it doesn't compile :

 real(wp), dimension(:,:), allocatable :: outpower
       1
Error: Symbol ‘wp’ at (1) has no IMPLICIT type
../results.f90:3660:25:

             allocate(outpower(points,ncol))
                     1
Error: Allocate-object at (1) is neither a data pointer nor an     allocatable variable
../results.f90:3676:16:

outpower(:,1) = exp(PK_data%matpower(:,1))
            1
Error: Unclassifiable statement at (1)
../results.f90:3679:20:

                 outpower(:,3) = exp(PK_data%vvpower(:,1))
                1
Error: Unclassifiable statement at (1)
compilation terminated due to -fmax-errors=4.
make[1]: *** [results.o] Error 1
make: *** [camb] Error 2

2) Also, I tried :

!Export files of total  matter power spectra in h^{-1} Mpc units, against k/h.
Type(MatterTransferData), intent(in) :: MTrans
Type(CAMBdata) :: State
character(LEN=Ini_max_string_len), intent(IN) :: FileNames(*)
character(LEN=name_tag_len) :: columns(3)
integer itf, i, unit
integer points
! Added : way of declaring double precision
double precision, dimension(:,:), allocatable :: outpower

but same thing, no compilation succeeded

call Transfer_GetMatterPowerS(State, MTrans, outpower(1,1), itf,  minkh,dlnkh, points)
                                                                                                      1
Error: Type mismatch in argument ‘outpower’ at (1); passed REAL(8) to    REAL(4)
make[1]: *** [results.o] Error 1
make: *** [camb] Error 2

UPDATE 1:

with -fmax-errors=1, I get the following :

 call Transfer_GetMatterPowerS(State, MTrans, outpower(1,1), itf,  minkh,dlnkh, points)
                                                                                                      1
Error: Type mismatch in argument ‘outpower’ at (1); passed REAL(8) to REAL(4)
compilation terminated due to -fmax-errors=1.

Except the solution given by @Steve with compilation option -freal-4-real-8, isn't really there another solution that I could include directly into code, i.e the section that I have given ?

UPDATE 2: here below the 3 relevant subroutines Transfer_GetMatterPowerS , Transfer_GetMatterPowerData and Transfer_SaveMatterPower that produces the error when trying to get double precision :

subroutine Transfer_GetMatterPowerS(State, MTrans, outpower, itf, minkh, dlnkh, npoints, var1, var2)
class(CAMBdata) :: state
Type(MatterTransferData), intent(in) :: MTrans
integer, intent(in) :: itf, npoints
integer, intent(in), optional :: var1, var2
real, intent(out) :: outpower(*)
real, intent(in) :: minkh, dlnkh
real(dl) :: outpowerd(npoints)
real(dl):: minkhd, dlnkhd

minkhd = minkh; dlnkhd = dlnkh
call Transfer_GetMatterPowerD(State, MTrans,  outpowerd, itf, minkhd, dlnkhd, npoints,var1, var2)
outpower(1:npoints) = outpowerd(1:npoints)

end subroutine Transfer_GetMatterPowerS

subroutine Transfer_GetMatterPowerData(State, MTrans, PK_data, itf_only, var1, var2)
!Does *NOT* include non-linear corrections
!Get total matter power spectrum in units of (h Mpc^{-1})^3 ready for interpolation.
!Here there definition is < Delta^2(x) > = 1/(2 pi)^3 int d^3k P_k(k)
!We are assuming that Cls are generated so any baryonic wiggles are well sampled and that matter power
!spectrum is generated to beyond the CMB k_max
class(CAMBdata) :: State
Type(MatterTransferData), intent(in) :: MTrans
Type(MatterPowerData) :: PK_data
integer, intent(in), optional :: itf_only
integer, intent(in), optional :: var1, var2
double precision :: h, kh, k, power
integer :: ik, nz, itf, itf_start, itf_end, s1, s2

s1 = PresentDefault (transfer_power_var, var1)
s2 = PresentDefault (transfer_power_var, var2)

if (present(itf_only)) then
    itf_start=itf_only
    itf_end = itf_only
    nz = 1
else
    itf_start=1
    nz= size(MTrans%TransferData,3)
    itf_end = nz
end if
PK_data%num_k = MTrans%num_q_trans
PK_Data%num_z = nz

allocate(PK_data%matpower(PK_data%num_k,nz))
allocate(PK_data%ddmat(PK_data%num_k,nz))
allocate(PK_data%nonlin_ratio(PK_data%num_k,nz))
allocate(PK_data%log_kh(PK_data%num_k))
allocate(PK_data%redshifts(nz))
PK_data%redshifts = State%Transfer_Redshifts(itf_start:itf_end)

h = State%CP%H0/100

do ik=1,MTrans%num_q_trans
    kh = MTrans%TransferData(Transfer_kh,ik,1)
    k = kh*h
    PK_data%log_kh(ik) = log(kh)
    power = State%CP%InitPower%ScalarPower(k)
    if (global_error_flag/=0) then
        call MatterPowerdata_Free(PK_data)
        return
    end if
    do itf = 1, nz
        PK_data%matpower(ik,itf) = &
            log(MTrans%TransferData(s1,ik,itf_start+itf-1)*&
            MTrans%TransferData(s2,ik,itf_start+itf-1)*k &
            *const_pi*const_twopi*h**3*power)
    end do
end do

call MatterPowerdata_getsplines(PK_data)

end subroutine Transfer_GetMatterPowerData

subroutine Transfer_SaveMatterPower(MTrans, State,FileNames, all21cm)
use constants

!Export files of total  matter power spectra in h^{-1} Mpc units, against k/h.
Type(MatterTransferData), intent(in) :: MTrans
Type(CAMBdata) :: State
character(LEN=Ini_max_string_len), intent(IN) :: FileNames(*)
character(LEN=name_tag_len) :: columns(3)
integer itf, i, unit
integer points
! Added : way of declaring double precision
!integer, parameter :: wp = selected_real_kind(15,307)
!real(wp), dimension(:,:), allocatable :: outpower
double precision, dimension(:,:), allocatable :: outpower
real minkh,dlnkh
Type(MatterPowerData) :: PK_data
integer ncol
logical, intent(in), optional :: all21cm
logical all21
!JD 08/13 Changes in here to PK arrays and variables
integer itf_PK

all21 = DefaultFalse(all21cm)
if (all21) then
    ncol = 3
else
    ncol = 1
end if

do itf=1, State%CP%Transfer%PK_num_redshifts
    if (FileNames(itf) /= '') then
        if (.not. transfer_interp_matterpower ) then
            itf_PK = State%PK_redshifts_index(itf)

            points = MTrans%num_q_trans
            allocate(outpower(points,ncol))

            !Sources
            if (all21) then
                call Transfer_Get21cmPowerData(MTrans, State, PK_data, itf_PK)
            else
                call Transfer_GetMatterPowerData(State, MTrans, PK_data, itf_PK)
                !JD 08/13 for nonlinear lensing of CMB + LSS compatibility
                !Changed (CP%NonLinear/=NonLinear_None) to CP%NonLinear/=NonLinear_none .and. CP%NonLinear/=NonLinear_Lens)
                if(State%CP%NonLinear/=NonLinear_none .and. State%CP%NonLinear/=NonLinear_Lens) then
                    call State%CP%NonLinearModel%GetNonLinRatios(State, PK_data)
                    PK_data%matpower = PK_data%matpower +  2*log(PK_data%nonlin_ratio)
                    call MatterPowerdata_getsplines(PK_data)
                end if
            end if

            outpower(:,1) = exp(PK_data%matpower(:,1))
            !Sources
            if (all21) then
                outpower(:,3) = exp(PK_data%vvpower(:,1))
                outpower(:,2) = exp(PK_data%vdpower(:,1))

                outpower(:,1) = outpower(:,1)/1d10*const_pi*const_twopi/MTrans%TransferData(Transfer_kh,:,1)**3
                outpower(:,2) = outpower(:,2)/1d10*const_pi*const_twopi/MTrans%TransferData(Transfer_kh,:,1)**3
                outpower(:,3) = outpower(:,3)/1d10*const_pi*const_twopi/MTrans%TransferData(Transfer_kh,:,1)**3
            end if

            call MatterPowerdata_Free(PK_Data)
            columns = ['P   ', 'P_vd','P_vv']
            unit = open_file_header(FileNames(itf), 'k/h', columns(:ncol), 15)
            do i=1,points
                write (unit, '(*(E15.6))') MTrans%TransferData(Transfer_kh,i,1),outpower(i,:)
            end do
            close(unit)
        else
            if (all21) stop 'Transfer_SaveMatterPower: if output all assume not interpolated'
            minkh = 1e-4
            dlnkh = 0.02
            points = log(MTrans%TransferData(Transfer_kh,MTrans%num_q_trans,itf)/minkh)/dlnkh+1
            !             dlnkh = log(MTrans%TransferData(Transfer_kh,MTrans%num_q_trans,itf)/minkh)/(points-0.999)
            allocate(outpower(points,1))
            call Transfer_GetMatterPowerS(State, MTrans, outpower(1,1), itf,  minkh,dlnkh, points)

            columns(1) = 'P'
            unit = open_file_header(FileNames(itf), 'k/h', columns(:1), 15)

            do i=1,points
                write (unit, '(*(E15.6))') minkh*exp((i-1)*dlnkh),outpower(i,1)
            end do
            close(unit)
        end if

        deallocate(outpower)
    end if
end do

end subroutine Transfer_SaveMatterPower

@VladimirF I can't make work the code here (sources are too big, that's why I put the relevant subsroutine), I can only show you the data generated and what the data I would like to get (see begin of post added). — youpilat13, Jun 04 '19 at 21:49
All the errors after the first one are run-on errors. Compile with -fmax-errors=1, and fix the error. You can try to skip a proper porting to double precision and use the -freal-4-real-8 option. — Steve, Jun 04 '19 at 22:55
You didn't put all the *relevant*.code, that's why I askwed. And I did not sask for all your ask, I asked to ake a [mcve], please do read the link. — Vladimir F, Jun 05 '19 at 05:50
I saw your UPDATE1. It appears you've changed the precision of one of the actual arguments to the function but failed to change the precision of the corresponding dummy argument. That's what the error message is telling you. I don't understand your comment about -freal-4-real-8 (unless your version of the compiler is ancient). — Steve, Jun 05 '19 at 17:54
@Steve . I would like to find a solution that avoids to use compilation flags to get double precision : I would like to implement this double precision directly in the code : from you experience, have you got any suggestions or tracks to give me to perform this modification ? — youpilat13, Jun 05 '19 at 19:48
I commend you not wanting to use an option. In fact, I'm fairly adamant about doing a proper porting job over a crutch option. For the job at hand, use -fmax-errors=1 and address each error as it occurs. In your Update 1, the error points to the `Transfer_GetMatterPowerS` subroutine and in particular the dummy assocated with the actual argument `outpower` has a type mismatch. You need to edit the subroutine and promote that dummy argument to double precision. Easiest way to do this is with a `types` module that is included where needed and changing `REAL` to `REAL(dp)` everywhere. — Steve, Jun 05 '19 at 20:55
@Steve thank you for your advice, I got it by replacing all `REAL` by `REAL(dp)`. Regards — youpilat13, Jun 06 '19 at 08:03

Can't replace Fortran real variables by double precision variables or more precision

0 Answers0