Can Fortran 77 Code Be Used to Debug Python Code for Solving ODEs Using Radau5?

Vick · May 14, 2023

Mark44 said:

In one of the subroutines there is this monstrosity

That is precisely the formulation: system of 66 1st Order Ordinary Differential Equations (ODEs) which I need, to solve the EMEP Problem. But as you can see the 66 ODEs have several variables (RC, DJ, Y, etc) for which I need to make sure I have the correct values.

Vick · May 14, 2023

PeterDonis said:

There is a website address in the code

It refers to the same web address I took the code from: the name of the web address being changed since the code was made!

Vick · May 14, 2023

PeterDonis said:

This appears to be just one problem in a test set. That might be why the FORTRAN code is not a self-contained program: it is intended to be run by some kind of test-running program that is somewhere else on the website you linked to. Again, the people who run the website would be much better able to clarify things like that than we are.

Yes I have my Python program solver to solve for stiff ODEs. But I need a formulation of the problem in Python. The formulation's source happens to be that Fortran code!!

Tom.G · May 14, 2023

The code you are working with is about 17 years old. The links at the top of the listing are dead, those pages no longer exist.

Since you obtained the code from an archive site, I looked there and found:

https://archimede.uniba.it/~testset/report/emep.pdf

Look on page 2 of that PDF file for "References." That will at least give you leads to follow in your search for more support information.

Good Luck!
Tom

FactChecker · May 14, 2023

The code in your post #35 should be a big help if someone wants to run the FORTRAN to mimic your Python. Print statements can be inserted to match the prints in your code.

PS. I actually found it much easier to download and look at your code in the way you included it in your first post. The copy and paste from post #35 was a harder process.

Vick · May 14, 2023

FactChecker said:

The copy and paste from post #35 was a harder process.

Yep!

Steve4Physics · May 14, 2023

Vick said:

TL;DR Summary: Debug Fortran Code

I have a piece of Fortran code but I'm not Fortran literate. I'm trying to translate it into Python. I have already made it, but I need someone who can run the Fortran code to compare the values of the variables.

Hi. I will chip in if I may. Many years (err, decades) ago I did a lot of scientific/technical programming in Fortran, so thought I’d take a look at this thread.

As already noted, the Fortran code (supplied in Post #1) is essentially only a collection of subroutines.

Each subroutine can only be called if you know appropriate values (or variable names) for its arguments (the input parameters). I can’t see how this would be possible without a great deal more information.

Of course, you’d also need to know in what order to call the subroutines. And if some sort of iterative loop were involved in the main program, you’d need to know what specific subroutines were in the loop.

Unless you can get your hands on an old copy of the main program, my opinion is that this is a lost cause.

[Minor edit.]

Vick · May 14, 2023

Steve4Physics said:

get your hands on an old copy of the main program

The main program is probably a Fortran stiff ODE solver. But my question is about formulating this code in Python because my "main program" is a Python stiff ODE solver. From my Python version, doesn't it show what are the values and variables that are used in order to call up some Fortran functions? I actually do not require to run the Fortran subroutines. I just need to know the values that the variables do take for a specified TIME.

TeethWhitener · May 14, 2023

I’ve only looked at the Fortran code very briefly, but this might be a partial answer.

Since the emepcf subroutine doesn’t seem to have any dependencies, you should be able to basically run it as a standalone program. Just delete the line that begins the subroutine and replace it with “program emepcf”. Then you can initialize the time where it says “double precision time” etc. Just put “double precision time=331000” or whatever the number you wanted to use was. You’ll also need to include a line at the end that outputs the DJ and RC variables. Any online fortran resource will tell you how to do that. That should give you DJ and RC for a certain time which you can compare to the output of your python code.

Getting the other variables you want from the other subroutines that call this one is a bit more involved, but if I were trying to debug this (or translate it to python, as the case may be), I’d probably follow the same principles (turn the subroutines into the main program, etc.)

Vick · May 14, 2023

TeethWhitener said:

Just delete the line that begins the subroutine and replace it with “program emepcf”. Then you can initialize the time where it says “double precision time” etc.

I don't have a Fortran compiler as I'm not Fortran literate. That's why I was asking help from someone who can do what you suggested and provide me with the values of the variables. I have already made a Python version of that code as best as I could understand from the math parts, and I wanted to compare values of variables.

Vick · May 14, 2023

TeethWhitener said:

Getting the other variables you want from the other subroutines that call this one is a bit more involved

You only need to get the values for the T, Y, DY, DJ and RC. You don't need to find for JAC (that's the Jacobian matrix which my solver computes automatically) and just below the JAC in the code are the results from their Fortran solver using Radau5 with different parameters. We don't need these either.

FactChecker · May 14, 2023

Steve4Physics said:

Hi. I will chip in if I may. Many years (err, decades) ago I did a lot of scientific/technical programming in Fortran, so thought I’d take a look at this thread.

As already noted, the Fortran code (supplied in Post #1) is essentially only a collection of subroutines.

Each subroutine can only be called if you know appropriate values (or variable names) for its arguments (the input parameters). I can’t see how this would be possible without a great deal more information.

Of course, you’d also need to know in what order to call the subroutines. And if some sort of iterative loop were involved in the main program, you’d need to know what specific subroutines were in the loop.

Unless you can get your hands on an old copy of the main program, my opinion is that this is a lost cause.

[Minor edit.]

I think the first objective is to verify that his translation of the FORTRAN to Python is correct. For that, arranging the FORTRAN code to do exactly what the Python code does would allow him to compare the (possibly artificial) calculation results. That is how the code in post #35 can be helpful. It shows what the Python program does and what FORTRAN code matches it. With some work, it should be possible to make a FORTRAN program that mimics the Python program. Write statements can be added to the FORTRAN that correspond to the prints in the Python program.

PeterDonis · May 14, 2023

Vick said:

I actually do not require to run the Fortran subroutines. I just need to know the values that the variables do take for a specified TIME.

I don't see how you would be able to do that without running the FORTRAN code and keeping track of its variable values.

TeethWhitener · May 14, 2023

Vick said:

I don't have a Fortran compiler

gfortran is free and easy to use. It’s part of the gcc system if you have that installed already.

Vick · May 14, 2023

PeterDonis said:

I don't see how you would be able to do that without running the FORTRAN code and keeping track of its variable values.

Python:

def test(a):
    c = a*4
    return c# without running test
c = a*4
print(c)

Vick · May 14, 2023

TeethWhitener said:

gfortran is free and easy to use. It’s part of the gcc system if you have that installed already.

I am not a Fortran user. You are asking me to install a program that I don't know how to use. Whereas you already know it and have it.

FactChecker · May 14, 2023

PeterDonis said:

I don't see how you would be able to do that without running the FORTRAN code and keeping track of its variable values.

I think he just needs data to verify/debug his conversion to Python code. He is looking for someone to run the FORTRAN and give him the results. In reality, the FORTRAN program needs to be developed to mimic his Python order of execution and outputs. I assume that he is knowledgeable enough in the subject matter to know the basic code organization that he needs for his use of the Python program. It's all the details of the calculation that he wants to verify. In post #35 he has inserted the FORTRAN code in comments that match the appropriate parts of his Python code.

PeterDonis · May 14, 2023

FactChecker said:

I think he just needs data to verify/debug his conversion to Python code. He is looking for someone to run the FORTRAN and give him the results.

Yes: someone to run the FORTRAN. In other words, the data he needs can only be obtained by running the FORTRAN code. So when he says, as in the post that I quoted, that "I actually do not require to run the Fortran subroutines", I don't see how that can be true: as you agree, someone needs to run the FORTRAN code to give him the data he needs.

TeethWhitener · May 14, 2023

This is turning into r/choosingbeggars

PeterDonis · May 14, 2023

Vick said:

I am not a Fortran user. You are asking me to install a program that I don't know how to use. Whereas you already know it and have it.

Fortran is not a "program". It's a programming language. The fact that someone knows that programming language does not mean that person knows how to write a FORTRAN program that will run the FORTRAN subroutines in the code you posted in the proper way to give the proper outputs that you are asking for.

Vanadium 50 · May 14, 2023

Vick said:

You are asking me to install a program that I don't know how to use.

And you are asking us - or to use your own words, "expecting us" - to debug non-working FORTRAN code for you. (At least two problems - no main program, and doesn't print the variables you want printed)

Who is being more reasonable?

berkeman · May 14, 2023

Thread closed temporarily for Moderation...

berkeman · May 15, 2023

After a Mentor discussion, this thread will remain closed. Asking us to do all of that work for the OP is too big of an ask. Thanks to all who gave the OP ideas for how to do this themself.

PeterDonis · May 16, 2023

Moderator's note: Thread reopened to allow @Vick to post an update.

Vick · May 16, 2023

With the suggestions from posts #2 and #44, I modified the code myself as per #44 and compiled it as per #2: the results was as I expected, viz. I did not quite fully understand the Fortran syntax even though, most of it is math. I had two errors when compared to my Python formulation.
I'd like to thank @TeethWhitener from post #44 for the suggestion and @Filip Larsen from post #2 for his suggestion. And also I'd like to thank @FactChecker for trying to understand my question from my perspective.
Below the modified Fortran code to print values of variables:

Fortran:

Program Hello
DOUBLE PRECISION Y(66), DY(66)
DOUBLE PRECISION HMIX
DOUBLE PRECISION M, O2, XN2, RPATH3, RPATH4, DELTA
DOUBLE PRECISION S1, S2, S3, S4, S5, S6, S7, S8, S9, S10, S11
DOUBLE PRECISION EMIS1, EMIS2, EMIS3, EMIS4, EMIS5, EMIS6,EMIS7, EMIS8, EMIS9, EMIS10, EMIS11, EMIS12, EMIS13
DOUBLE PRECISION FRAC6, FRAC7, FRAC8, FRAC9, FRAC10, FRAC13
DOUBLE PRECISION VEKT6, VEKT7, VEKT8, VEKT9, VEKT10, VEKT13
DOUBLE PRECISION VMHC, EMNOX, EMHC, FACISO, FACHC, FACNOX
DOUBLE PRECISION RC(266), DJ(16), H2O
DOUBLE PRECISION A(16), B(16), SEC, T
INTEGER ITIMEH, I24HRS, I
DOUBLE PRECISION TIME, TIMEH, TIMEOD, PI, XLHA, FI, DEKL, XQ, RH, XZDO I = 1, 13
Y(I)=1.D7
end DO

DO I = 14, 66
Y(I)=100.D0
end DO
Y(1) = 1.0D9
Y(2) = 5.0D9
Y(3) = 5.0D9
Y(4) =3.8D12
Y(5) =3.5D13
Y(14)=5.D11
Y(38)=1.D-3HMIX=1.2D5
FRAC6=0.07689D0
FRAC7=0.41444D0
FRAC8=0.03642D0
FRAC9=0.03827D0
FRAC10=0.24537D0
FRAC13=0.13957D0
VEKT6=30.D0
VEKT7=58.D0
VEKT8=28.D0
VEKT9=42.D0
VEKT10=106.D0
VEKT13=46.D0
VMHC=1.D0/(FRAC6/VEKT6+FRAC7/VEKT7+FRAC8/VEKT8+FRAC9/VEKT9+FRAC10/VEKT10+FRAC13/VEKT13)

EMNOX=2.5D11
EMHC=2.5D11

FACISO=1.0D0
FACHC=1.0D0
FACNOX=1.0D0

EMIS1 = EMNOX * FACNOX
EMIS2 = 0.D0
EMIS3 = EMNOX * FACNOX
EMIS4 = EMHC*10.D0 * FACHC
EMIS5 = 0.D0
EMIS6 = EMHC*FRAC6/VEKT6*VMHC  * FACHC
EMIS7 = EMHC*FRAC7/VEKT7*VMHC  * FACHC
EMIS8 = EMHC*FRAC8/VEKT8*VMHC  * FACHC
EMIS9 = EMHC*FRAC9/VEKT9*VMHC  * FACHC
EMIS10= EMHC*FRAC10/VEKT10*VMHC* FACHC
EMIS11= 0.5D0 * FACISO * EMHC
EMIS12= 0.D0
EMIS13= EMHC*FRAC13/VEKT13*VMHC* FACHC

M = 2.55D19
O2 = 5.2D18
XN2= 1.99D19

RPATH3 = 0.65D0
RPATH4 = 0.35D0A = (/1.23D-3, 2.00D-4, 1.45D-2, 2.20D-5, 3.00D-6, 5.40D-5, &
6.65D-5, 1.35D-5, 2.43D-5, 5.40D-4, 2.16D-4, 5.40D-5, 3.53D-2, 8.94D-2, 3.32D-5, 2.27D-5/)

B = (/0.60D0, 1.40D0, 0.40D0, 0.75D0, 1.25D0, 0.79D0, 0.60D0, &
0.94D0, 0.88D0, 0.79D0, 0.79D0, 0.79D0, 0.081D0, 0.059D0, 0.57D0, 0.62D0/)

TIME = 325025
TIMEH=TIME/3600.D0
ITIMEH=int(TIMEH)
I24HRS=ITIMEH/24+1
TIMEOD=TIMEH-(I24HRS-1)*24.D0

PI=4.D0*ATAN(1.0D0)

XLHA=(1.D0+TIMEOD*3600.D0/4.32D4)*PI

FI=50.D0*PI/180.D0
DEKL=23.5D0*PI/180.D0
SEC=1.D0/(COS(XLHA)*COS(FI)*COS(DEKL)+SIN(FI)*SIN(DEKL))

T=298.D0

XQ=-7.93D-5*TIMEOD*3600.D0+2.43D0
RH=23.D0*SIN(XQ)+66.5D0

XZ=(597.3D0-0.57D0*(T-273.16D0))*18.D0/1.986D0*(1.D0/T-1.D0/273.16D0)
H2O=6.1078D0*EXP(-XZ)*10.D0/(1.38D-16*T)*RH

IF (TIMEOD .LT. 4.0D0 .OR. TIMEOD .GE. 20.D0) THEN

         DO 100 I = 1,16
            DJ(I)=1.D-40
 100     CONTINUE
      ELSE

         DO 110 I = 1,16
            DJ(I)=A(I) * EXP( -B(I) * SEC )
            IF( DJ(I) .LT. 0.0D0 ) STOP 'DJ'
 110     CONTINUE
      ENDIF

     

DO 120 I = 1,266
    RC(I)=0.D0
 120  CONTINUE

    rc(1)  =5.7d-34*(t/300.0D0)**(-2.8D0) * m * o2

      rc(5)  =9.6d-32*(t/300.0D0)**(-1.6D0) * m

      rc(7)  =2.0d-11*exp(100.D0/t) * m

      rc(11) =1.8d-12*exp(-1370.D0/t)
      rc(12) =1.2d-13*exp(-2450.D0/t)
      rc(13) =1.9d-12*exp(-1000.D0/t)
      rc(14) =1.4d-14*exp(-600.D0/t)
      rc(15) =1.8d-11*exp(+110.D0/t)
      rc(17) =3.7d-12*exp(240.0D0/t)

      rc(19) =7.2d-14*exp(-1414.D0/t)

      rc(29) =7.1d14*exp(-11080.D0/t)

      rc(30) =4.8d-11*exp(+250.D0/t)

      rc(31) =2.9d-12*exp(-160.D0/t)

      rc(33) =7.7d-12*exp(-2100.D0/t)

      rc(35) =1.0d-14*exp(785.0D0/t)

      rc(36) =2.3d-13*exp(600.D0/t)
      rc(36) = rc(36) + m * 1.7d-33*exp(1000.D0/t)
      rc(36) = rc(36) *(1.D0 + 1.4d-21 * h2o *exp(2200.D0/t))      rc(59) =3.9d-12*exp(-1885.D0/t)

      rc(60) =4.2d-12*exp(180.D0/t)

      rc(61) =5.5d-14*exp(365.D0/t)

      rc(62) =5.5d-14*exp(365.D0/t)

      rc(63) =3.3d-12*exp(-380.D0/t)

      rc(65) =3.8d-13*exp(780.D0/t)

      rc(67) =1.0d-12*exp(190.D0/t)

      rc(68) =1.9d-12*exp(190.D0/t)

      rc(71) =7.8d-12*exp(-1020.D0/t)

      rc(74) =6.5d-13*exp(650.D0/t)

      rc(75) =5.6d-12*exp(310.D0/t)

      rc(76) = 5.8d-12*exp(190.D0/t)

      rc(78) =1.34d16*exp(-13330.D0/t)

      rc(88) =1.3d-13*exp(1040.D0/t)

      rc(89) =3.0d-13*exp(1040.D0/t)

      rc(94) =2.8d-12*exp(530.D0/t)

      rc(81) =1.64d-11*exp(-559.D0/t)      rc(83) =rc(60)
      rc(105)=rc(60)

      rc(110)=rc(60)

      rc(162)=rc(60)

      rc(164)=rc(60)      rc(109)=1.66d-12*exp(474.D0/t)

      rc(112)=1.2d-14*exp(-2630.D0/t)

      rc(123)=6.5d-15*exp(-1880.D0/t)

      rc(126) = rc(60)

      rc(220) = rc(60)

      rc(236) = rc(60)

      rc(150) = 12.3d-15*exp(-2013.D0/t)

      rc(151) = 2.54d-11*exp(410.D0/t)

      rc(152) = rc(60)

      rc(153) = 4.13d-12*exp(452.D0/t)

      rc(154) = rc(60)

      rc(158) = 1.86d-11*exp(175.D0/t)

      rc(159) = rc(60)

      rc(160) = 4.32d-15*exp(-2016.D0/t)

      rc(43)=5.d-6
      if(rh.gt.0.90D0) rc(43)=1.d-4
      rc(44) =rc(43)
      rc(45) =rc(43)      rc(8)  =2.2d-10

      rc(20) =1.4d-12

      rc(21) =1.4d-11

      rc(26) =4.1d-16

      rc(39) =1.35d-12

      rc(40) =4.0d-17

      rc(64) =3.2d-12

      rc(66) =9.6d-12

      rc(69) =5.8d-16
      rc(70) =2.4d-13
      rc(72) =8.9d-12      rc(77) =1.0d-11

      rc(79) =2.0d-11

      rc(80) =1.1d-11

      rc(85) =1.0d-11

      rc(86) =1.15d-12

      rc(87)= 4.8d-12

      rc(90) =8.0d-13

      rc(125)=2.86d-11

      rc(146) = 3.2d-11

      rc(147) = 2.0d-11

      rc(148) = 2.2d-11

      rc(149) = 3.7d-11

      rc(155) =rc(85)

      rc(156) =2.0d-11

      rc(157) =8.0d-18

      rc(161) =3.35d-11

      rc(163) =7.8d-13

      rc(219)=2.0d-11

      rc(221)=1.1d-11

      rc(222)=1.7d-11

      rc(223)= 2.4d-11

      rc(234)=1.37d-11

      rc(235)= 1.7d-11

          delta=1.0D0
          if(timeod.ge.20.D0.or.timeod.lt.4.D0) delta=0.25D0

          rc(46) =2.0D0 * delta /hmix
          rc(47) =0.5D0 * delta /hmix
          rc(48) =0.2D0 * delta /hmix
          rc(49) =0.5D0 * delta /hmix
          rc(50) =0.2D0 * delta /hmix
          rc(51) =0.1D0 * delta /hmix

          rc(52) = 0.5D0 *delta /hmix

          rc(53) = 0.3D0 *delta /hmix
         
         
DY(1) = DJ(3)*Y(2)+DJ(13)*Y(39)+RC(19)*Y(2)*Y(39)+EMIS1/HMIX-&
(RC(5)*Y(36)+RC(11)*Y(14)+RC(17)*Y(15)+RC(72)*Y(23)+RC(79)*(Y(24)+Y(57))+&
RC(15)*Y(39)+RC(60)*(Y(19)+Y(26)+Y(27)+Y(29)+Y(31)+Y(33)+Y(35)+&
Y(43)+Y(45)+Y(59)+Y(61)+Y(60)))*Y(1)
DY(2) = Y(1)*(RC(5)*Y(36)+RC(11)*Y(14)+RC(17)*Y(15)+RC(72)*Y(23)+&
RC(79)*(Y(24)+Y(57))+0.2D1*RC(15)*Y(39))+RC(60)*Y(1)*(Y(19)+Y(26)+&
Y(27)+Y(29)+Y(31)+Y(33)+Y(35)+Y(59))+RC(60)*Y(1)*(0.86D0*Y(43)+0.19D1*&
Y(61)+0.11D1*Y(60)+0.95D0*Y(45))+DJ(14)*Y(39)+DJ(5)*Y(16)+DJ(15)*&
Y(40)+RC(29)*Y(40)+RC(78)*(Y(25)+Y(58))-(DJ(3)+RC(12)*Y(14)+RC(20)*&
Y(39)+RC(21)*Y(37)+RC(48)+RC(77)*(Y(24)+Y(57)))*Y(2)
DY(3) = EMIS3/HMIX-(RC(39)*Y(37)+RC(40)*Y(19)+RC(47))*Y(3)
DY(4) = EMIS4/HMIX+Y(37)*(RC(66)*Y(11)+2.D0*RC(221)*Y(32)+RC(222)*&
Y(30))+Y(14)*(0.44D0*RC(112)*Y(8)+0.4D0*RC(123)*Y(9)+0.5D-1*RC(160)*&
Y(44)+0.5D-1*RC(150)*Y(41))+Y(11)*(DJ(6)+DJ(7)+RC(69)*Y(39))+DJ(8)*&
Y(12)+DJ(11)*Y(30)+2.D0*DJ(7)*Y(32)-RC(70)*Y(37)*Y(4)
DY(5) = EMIS5/HMIX+0.7D-1*RC(123)*Y(14)*Y(9)-RC(59)*Y(37)*Y(5)
DY(6) = EMIS6/HMIX-RC(71)*Y(37)*Y(6)
DY(7) = EMIS7/HMIX-RC(81)*Y(37)*Y(7)
DY(8) = EMIS8/HMIX-(RC(109)*Y(37)+RC(112)*Y(14))*Y(8)
DY(9) = EMIS9/HMIX+0.7D-1*RC(150)*Y(14)*Y(41)-(RC(123)*Y(14)+RC(125)*Y(37))*Y(9)
DY(10) = EMIS10/HMIX-RC(234)*Y(37)*Y(10)
DY(11) = Y(19)*(RC(60)*Y(1)+(2.D0*RC(61)+RC(62))*Y(19)+RC(80)*Y(24)+&
RC(40)*Y(3))+Y(37)*(RC(63)*Y(46)+RC(67)*Y(22))+Y(1)*RC(60)*(2.D0*&
Y(29)+Y(31)+0.74D0*Y(43)+0.266D0*Y(45)+0.15D0*Y(60))+Y(14)*(0.5D0*&
RC(123)*Y(9)+RC(112)*Y(8)+0.7D0*RC(157)*Y(56)+0.8D0*RC(160)*Y(44)+&
0.8D0*RC(150)*Y(41))+2.D0*DJ(7)*Y(32)+DJ(16)*(Y(22)+0.156D1*Y(50)+&
Y(51))-(RC(66)*Y(37)+DJ(6)+DJ(7)+RC(69)*Y(39)+RC(53))*Y(11)
DY(12) = Y(1)*(RC(72)*Y(23)+RC(83)*Y(26)*RPATH4+RC(105)*Y(27)+RC(126)*&
Y(31)+0.95D0*RC(162)*Y(61)+0.684D0*RC(154)*Y(45))+Y(37)*(RC(64)*&
Y(20)+RC(76)*Y(28)+RC(76)*Y(50))+0.5D0*RC(123)*Y(14)*Y(9)+0.4D-1*&
RC(160)*Y(14)*Y(44)+DJ(16)*(Y(28)+0.22D0*Y(50)+0.35D0*Y(49)+Y(51)+&
Y(52))-(DJ(8)+RC(75)*Y(37)+RC(53))*Y(12)
DY(13) = RC(83)*Y(1)*Y(26)*RPATH3+(0.65D0*DJ(16)+RC(76)*Y(37))*Y(49)+&
RC(76)*Y(37)*Y(51)+(RC(159)*Y(1)+DJ(16))*Y(59)+0.95D0*RC(162)*Y(1)*&
Y(61)-(DJ(9)+RC(86)*Y(37)+RC(53))*Y(13)
DY(14) = RC(1)*Y(36)+RC(89)*Y(15)*Y(24)-(RC(11)*Y(1)+RC(12)*Y(2)+RC(13)*&
Y(37)+RC(14)*Y(15)+RC(49)+RC(112)*Y(8)+RC(123)*Y(9)+RC(157)*&
Y(56)+RC(160)*Y(44)+RC(150)*Y(41)+DJ(1)+DJ(2))*Y(14)
s1 = Y(37)*(RC(13)*Y(14)+RC(31)*Y(17)+RC(33)*Y(18)+RC(39)*Y(3)+RC(63)*&
Y(46)+RC(64)*Y(20)+RC(66)*Y(11)+RC(70)*Y(4)+RC(221)*Y(32))+Y(19)*&
(RC(40)*Y(3)+2.D0*RC(61)*Y(19)+0.5D0*RC(80)*Y(24))+DJ(11)*Y(30)+&
Y(1)*RC(60)*(Y(19)+Y(29)+Y(31)+Y(33)+Y(35)+0.95D0*Y(45)+Y(26)*RPATH3+&
0.78D0*Y(43)+Y(59)+0.5D-1*Y(61)+0.8D0*Y(60))+RC(72)*Y(1)*Y(23)+DJ(8)*Y(12)
DY(15) = s1+2.D0*DJ(6)*Y(11)+DJ(16)*(Y(22)+Y(28)+0.65D0*Y(49)+Y(50)+&
Y(51)+Y(48)+Y(53))+Y(39)*(RC(26)*Y(17)+RC(69)*Y(11))+Y(14)*(0.12D0*&
RC(112)*Y(8)+0.28D0*RC(123)*Y(9)+0.6D-1*RC(160)*Y(44))+0.6D-1*&
RC(150)*Y(14)*Y(41)-(RC(14)*Y(14)+RC(17)*Y(1)+RC(30)*Y(37)+2.D0*RC(36)*&
Y(15)+RC(65)*Y(19)+RC(74)*Y(23)+(RC(88)+RC(89))*Y(24)+RC(85)*&
(Y(26)+Y(29)+Y(31)+Y(27)+Y(57)+Y(45)+Y(61)+Y(59)+Y(33)+Y(35)+Y(43)+&
Y(60)))*Y(15)
DY(16) = RC(21)*Y(2)*Y(37)+Y(39)*(RC(26)*Y(17)+RC(69)*Y(11))-(RC(35)*Y(37)+DJ(5)+RC(45))*Y(16)
DY(17) = RC(36)*Y(15)**2-(RC(31)*Y(37)+DJ(4)+RC(43)+RC(26)*Y(39)+RC(47))*Y(17)
DY(18) = DJ(7)*Y(11)+Y(14)*(0.13D0*RC(112)*Y(8)+0.7D-1*RC(123)*Y(9))-RC(33)*Y(37)*Y(18)
DY(19) = Y(37)*(RC(59)*Y(5)+RC(68)*Y(22))+Y(24)*(RC(79)*Y(1)+2.D0*&
RC(94)*Y(24))+DJ(8)*Y(12)+DJ(16)*Y(47)+0.31D0*RC(123)*Y(14)*Y(9)-&
(RC(40)*Y(3)+RC(60)*Y(1)+2.D0*RC(61)*Y(19)+2.D0*RC(62)*Y(19)+RC(65)*&
Y(15)+0.5D0*RC(80)*Y(24))*Y(19)
DY(20) = EMIS13/HMIX-RC(64)*Y(37)*Y(20)
DY(21) = (RC(40)*Y(19)+RC(39)*Y(37))*Y(3)+0.5D-1*EMIS3/HMIX-RC(51)*Y(21)
DY(22) = RC(65)*Y(15)*Y(19)-(RC(43)+DJ(16)+(RC(67)+RC(68))*Y(37))*Y(22)
DY(23) = Y(37)*(RC(71)*Y(6)+RC(68)*Y(28))+0.35D0*DJ(16)*Y(49)+RC(83)*Y(1)*Y(26)*RPATH4+DJ(9)*Y(13)-(RC(72)*Y(1)+RC(74)*Y(15))*Y(23)
DY(24) = Y(37)*(RC(75)*Y(12)+RC(222)*Y(30)+RC(68)*Y(47))+RC(105)*&
Y(1)*Y(27)+RC(78)*Y(25)+DJ(11)*Y(30)+DJ(9)*Y(13)+DJ(16)*Y(52)+0.684D0*&
RC(154)*Y(1)*Y(45)-(RC(77)*Y(2)+RC(79)*Y(1)+RC(80)*Y(19)+2.D0*&
RC(94)*Y(24)+(RC(88)+RC(89))*Y(15))*Y(24)
DY(25) = RC(77)*Y(24)*Y(2)-(RC(50)+RC(78))*Y(25)
DY(26) = Y(37)*(RC(81)*Y(7)+RC(68)*Y(49))-(RC(83)*Y(1)+RC(85)*Y(15))*Y(26)
DY(27) = Y(37)*(RC(86)*Y(13)+RC(87)*Y(52))-(RC(105)*Y(1)+RC(85)*Y(15))*Y(27)
DY(28) = RC(74)*Y(15)*Y(23)-((RC(76)+RC(68))*Y(37)+DJ(16)+RC(52))*Y(28)
DY(29) = Y(37)*(RC(109)*Y(8)+RC(68)*Y(50))-(RC(110)*Y(1)+RC(85)*Y(15))*Y(29)
DY(30) = RC(236)*Y(1)*Y(33)+RC(220)*Y(1)*Y(35)+0.266D0*RC(154)*Y(1)*Y(45)+0.82D0*RC(160)*&
Y(14)*Y(44)+DJ(16)*(Y(48)+Y(53))-(DJ(11)+RC(222)*Y(37))*Y(30)
DY(31) = Y(37)*(RC(125)*Y(9)+RC(68)*Y(51))-(RC(126)*Y(1)+RC(85)*Y(15))*Y(31)
DY(32) = RC(220)*Y(1)*Y(35)+DJ(16)*Y(53)-(2.D0*DJ(7)+RC(221)*Y(37))*Y(32)
DY(33) = Y(37)*(RC(234)*Y(10)+RC(235)*Y(48))-(RC(236)*Y(1)+RC(85)*Y(15))*Y(33)
DY(34) = RC(236)*Y(1)*Y(33)+DJ(16)*Y(48)-RC(219)*Y(37)*Y(34)
DY(35) = Y(37)*(RC(219)*Y(34)+RC(223)*Y(53))-(RC(220)*Y(1)+RC(85)*Y(15))*Y(35)
DY(36) = DJ(1)*Y(14)+DJ(3)*Y(2)+DJ(14)*Y(39)+RC(7)*Y(38)+0.2D0*&
RC(160)*Y(14)*Y(44)+0.3D0*RC(150)*Y(14)*Y(41)-(RC(1)+RC(5)*Y(1))*Y(36)
s1 = 2.D0*RC(8)*H2O*Y(38)+Y(15)*(RC(14)*Y(14)+RC(17)*Y(1))+2.D0*DJ(4)*Y(17)+DJ(5)*Y(16)
s2 = s1+DJ(16)*(Y(22)+Y(28)+Y(47)+Y(49)+Y(50)+Y(52)+Y(48)+Y(53))
s3 = s2+Y(14)*(0.15D0*RC(123)*Y(9)+0.8D-1*RC(160)*Y(44))
s4 = s3
s6 = 0.55D0*RC(150)*Y(14)*Y(41)
s8 = -1
s11 = RC(222)*Y(30)+RC(75)*Y(12)+RC(81)*Y(7)+RC(87)*Y(52)+RC(86)*&
Y(13)+RC(235)*Y(48)+RC(109)*Y(8)+RC(125)*Y(9)+RC(234)*Y(10)+RC(223)*&
Y(53)+RC(219)*Y(34)+RC(31)*Y(17)+RC(21)*Y(2)+RC(148)*Y(62)+RC(64)*&
Y(20)+RC(39)*Y(3)+RC(71)*Y(6)
s10 = s11+RC(30)*Y(15)+RC(59)*Y(5)+RC(70)*Y(4)+RC(35)*Y(16)+RC(13)*&
Y(14)+RC(221)*Y(32)+RC(68)*(Y(22)+Y(28)+Y(47)+Y(50)+Y(51)+Y(49))+&
RC(66)*Y(11)+RC(151)*Y(41)+RC(153)*Y(44)+RC(63)*Y(46)+RC(33)*Y(18)+&
RC(158)*Y(54)+RC(146)*Y(63)+RC(149)*(Y(65)+Y(66))+RC(147)*Y(64)+RC(161)*Y(55)
s11 = Y(37)
s9 = s10*s11
s7 = s8*s9
s5 = s6+s7
DY(37) = s4+s5
DY(38) = DJ(2)*Y(14)-(RC(7)+RC(8)*H2O)*Y(38)
DY(39) = (RC(29)+DJ(15))*Y(40)+RC(12)*Y(14)*Y(2)+RC(35)*Y(37)*&
Y(16)-(RC(15)*Y(1)+RC(26)*Y(17)+RC(163)*Y(41)+RC(19)*Y(2)+RC(20)*Y(2)+DJ(13)+DJ(14)+RC(69)*Y(11))*Y(39)
DY(40) = RC(20)*Y(39)*Y(2)-(RC(29)+DJ(15)+RC(45))*Y(40)
DY(41) = EMIS11/HMIX-(RC(151)*Y(37)+RC(163)*Y(39)+RC(150)*Y(14))*Y(41)
DY(42) = RC(45)*Y(16)+2.D0*RC(44)*Y(40)-RC(51)*Y(42)
DY(43) = Y(37)*(RC(151)*Y(41)+RC(156)*Y(56))+0.12D0*RC(152)*Y(1)*Y(43)-(RC(152)*Y(1)+RC(155)*Y(15))*Y(43)
DY(44) = RC(60)*Y(1)*(0.42D0*Y(43)+0.5D-1*Y(60))+0.26D0*RC(150)*Y(14)*Y(41)-(RC(153)*Y(37)+RC(160)*Y(14))*Y(44)
DY(45) = RC(153)*Y(44)*Y(37)+RC(148)*Y(37)*Y(62)-(RC(154)*Y(1)+RC(85)*Y(15))*Y(45)
DY(46) = RC(62)*Y(19)**2-RC(63)*Y(37)*Y(46)
DY(47) = RC(88)*Y(15)*Y(24)-(RC(68)*Y(37)+DJ(16)+RC(52))*Y(47)
DY(48) = RC(85)*Y(15)*Y(33)-(RC(235)*Y(37)+DJ(16)+RC(52))*Y(48)
DY(49) = RC(85)*Y(15)*Y(26)-((RC(76)+RC(68))*Y(37)+DJ(16)+RC(52))*Y(49)
DY(50) = RC(85)*Y(15)*Y(29)-((RC(76)+RC(68))*Y(37)+DJ(16)+RC(52))*Y(50)
DY(51) = RC(85)*Y(15)*Y(31)-((RC(76)+RC(68))*Y(37)+DJ(16)+RC(52))*Y(51)
DY(52) = RC(85)*Y(15)*Y(27)-(RC(87)*Y(37)+DJ(16)+RC(52))*Y(52)
DY(53) = RC(85)*Y(15)*Y(35)-(RC(223)*Y(37)+DJ(16)+RC(52))*Y(53)
DY(54) = RC(60)*Y(1)*(0.32D0*Y(43)+0.1D0*Y(60))+0.67D0*RC(150)*Y(14)*Y(41)-RC(158)*Y(37)*Y(54)
DY(55) = RC(60)*Y(1)*(0.14D0*Y(43)+0.5D-1*Y(45)+0.85D0*Y(60))-RC(161)*Y(37)*Y(55)
DY(56) = RC(155)*Y(15)*Y(43)-(RC(156)*Y(37)+RC(157)*Y(14)+RC(52))*Y(56)
DY(57) = 0.5D0*RC(158)*Y(37)*Y(54)+RC(78)*Y(58)+RC(149)*Y(37)*Y(66)-(RC(77)*Y(2)+RC(79)*Y(1)+RC(85)*Y(15))*Y(57)
DY(58) = RC(77)*Y(57)*Y(2)-(RC(50)+RC(78))*Y(58)
DY(59) = RC(79)*Y(1)*Y(57)+RC(146)*Y(37)*Y(63)-(RC(159)*Y(1)+RC(85)*Y(15))*Y(59)
DY(60) = RC(163)*Y(39)*Y(41)+RC(147)*Y(37)*Y(64)-(RC(164)*Y(1)+RC(85)*Y(15))*Y(60)
DY(61) = RC(161)*Y(37)*Y(55)+RC(149)*Y(37)*Y(65)-(RC(162)*Y(1)+RC(85)*Y(15))*Y(61)
DY(62) = RC(85)*Y(15)*Y(45)-(RC(148)*Y(37)+RC(52))*Y(62)
DY(63) = RC(85)*Y(15)*Y(59)-(RC(146)*Y(37)+RC(52))*Y(63)
DY(64) = RC(85)*Y(15)*Y(60)-(RC(147)*Y(37)+RC(52))*Y(64)
DY(65) = RC(85)*Y(15)*Y(61)-(RC(149)*Y(37)+RC(52))*Y(65)
DY(66) = RC(85)*Y(15)*Y(57)-(RC(149)*Y(37)+RC(52))*Y(66)

Print *, DY
End Program Hello

Vick · May 16, 2023

PS:
Below the results from solving the 66 ODEs in Python using Radau5:

BvU · Feb 25, 2024

For what it's worth: I found the program that calls the routines in emep.f
The site has been moved to http://pitagora.dm.uniba.it/~testset/solvers/radau5.php
Pdf: radau5
The program is radau5d
It calls prob, init etc, and radau5
Other stuff in radaua and report

PDF http://archimede.dm.uniba.it/~testset/report/prologue.pdf has some documentation (page 16 onwards)

IV.3 Format of the problem codes
The eight subroutines that define the problem are called PROB, INIT, SETTOLERANCES, SETOUTPUT,FEVAL, JEVAL, MEVAL, and SOLUT. The following subsections describe the format of these subroutinesin full detail. An additional function PIDATE ...

Vick's emep.f from https://archimede.uniba.it/~testset/problems/emep.php

[edit]
@Vick: no fortran needed so far. I see that report.f has some utility to generate a
MATLAB and a SCILAB function to print the plots of the solution

just in case I can get a fortran compiler to run on this ancient PC: can you post a link to the .py source used for post #61 ?

##\ ##

Can Fortran 77 Code Be Used to Debug Python Code for Solving ODEs Using Radau5?

Similar threads

Hot Threads

Recent Insights