Replacement Patterns in FortranForm

Question

I have recently asked a question concerning string replacement rules for FortranForm. This question has been answered by Carl Woll - thanks for that!

Now as I was implementing this I was thinking of overcoming number definitions and simply replace real numbers with something like 1. --> 1._rkind. This replacement should work for all fractions, powers, etc.

I have found a post where this task has been tackled but unfortunately my Mathematica skills are very basic and I do not really understand what is going on there. Anyway, I have tried and this is what I was coming up so far:

G = -4/3 pi (-RAi + RAip1) (RAi^2 + RAi RAip1 + RAip1^2) rhoAi;
testReplace[someEq_] := Module[{trep},
   Clear[fort, fortranreal];
   fort[x_Real] := fortranreal[x];
   fort[x_] := x;
   trep = StringReplace[ToString@MapAll[fort, FortranForm[#]], Shortest["fortranreal(" ~~ s : __ ~~ ")"] :> s <> "_rkind"] & /@ someEq
];

This leaves me with the (unfortunately not quite useful) output:

(-4*pi*(-RAi + RAip1)*(RAi**2 + RAi*RAip1 + RAip1**2)*rhoAi)/3.

because obviously 4 is not an integer and would cause problems with the Fortran compiler. Next thing I've tried was replacing someEq with N[someEq] in the Module output which then leaves me with:

-1.3333333333333333_rkind pi RAi**2 + RAi*RAip1 + RAip1**2 RAip1 + RAi*-1._rkind rhoAi

where all the brackets are stripped away and e.g. -RAi gets replaced with RAi*-1._rkind.

What I am trying to achieve is evaluating all fractions and then attach _rkind to all real values. Can anyone tell me what am I doing wrong?

Answer 1

After running a few tests and digging deeper into Fortran standards (thanks @george2079 for pushing this into the right direction), I can finally post an aswer.

The method above is right but will cause problems when adding negative reals. As soon as we plug in e.g. 4./3. it will get converted to a real (1.33) and because the constant literal of 1.33 != 1.33_rkind the method attaches the _rkind to the number. Unfortunately Mathematica's Times[] will shuffle the number to the very end of the output. For positive values this does not matter but in case of negative values this will cause problems as we end up with e.g.:

G = -4./3 x;
testReplace[G]

Out:

x*-1.3333333333333_rkind >>> read update!

In order to fix this I've simply stopped Times[] by doing this with ClearAttributes[Times, Orderless]. So for completion the whole thing now looks like this:

Unprotect[Times];
ClearAttributes[Times, Orderless];
Protect[Times];

testReplace[someEq_] := Module[{trep}, Clear[fort, fortranreal];
  fort[x_Real] := fortranreal[x];
  fort[x_] := x;
  trep = StringReplace[ToString@MapAll[fort, FortranForm[#]],
    Shortest["fortranreal(" ~~ s : __ ~~ ")"] :> s <> "_rkind"] & /@ {someEq}[[1]]
];

G = -4./3 pi (-RAi + RAip1) (RAi^2 + RAi RAip1 + RAip1^2) rhoAi;
testReplace[G]

Out:

-1.333333333333333_rkind*pi*(-RAi + RAip1)*(RAi**2 + RAi*RAip1 + RAip1**2)*rhoAi

UPDATE

Actually you there is no need of fixing the reordering here as a*-b, a+-b and all other variations are allowed (thanks @george for the comment!) under Fortran (ifort and gfortran), it's just a bit better for readability.