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

In summary, the person is seeking help with translating a piece of Fortran code into Python, but they are not familiar with Fortran and need someone to run the code and provide the values of the variables for comparison. They mention that the code is a formulation for solving atmospheric emissions and that they want someone to run it on their compiler and provide the results. However, the other person in the conversation suggests that it may be difficult to find someone who is familiar with both FORTRAN and this specific code, and suggests looking for a person at the source of the code.
  • #36
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.
 
Technology news on Phys.org
  • #37
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!
 
  • #38
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!!
 
  • #39
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
 
  • #40
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.
 
  • Like
Likes Tom.G
  • #41
FactChecker said:
The copy and paste from post #35 was a harder process.
Yep!
 
  • #42
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.]
 
Last edited:
  • Like
Likes Vanadium 50
  • #43
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.
 
  • #44
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.)
 
  • Like
Likes Vick
  • #45
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.
 
  • #46
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.
 
  • #47
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.
 
  • Like
Likes Vick
  • #48
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.
 
  • Like
Likes Steve4Physics and TeethWhitener
  • #49
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.
 
  • #50
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)
 
  • #51
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.
 
  • #52
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.
 
  • Like
Likes Vick
  • #53
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.
 
  • #54
This is turning into r/choosingbeggars
 
  • #55
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.
 
  • #56
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?
 
  • #57
Thread closed temporarily for Moderation...
 
  • Like
Likes Tom.G
  • #58
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.
 
  • Like
Likes FactChecker
  • #60
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
 
  • Like
Likes Filip Larsen, pbuk and FactChecker
  • #61
PS:
Below the results from solving the 66 ODEs in Python using Radau5:
Figure_1.png
 
  • Like
Likes Tom.G, pbuk and FactChecker
  • #62
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 ?

##\ ##
 
Last edited:
  • #63
BvU said:
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)


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 ?

##\ ##
The testset emep.f is Fortran coding. As I am using Python, I needed a translation of the emep.f into Python. I succeeded and post 61 is the result of my own code using my own Radau5 code. Thanks
 
  • Like
Likes BvU

Similar threads

Fortran Porting Python Code to Fortran for Parallel Computing
    • Programming and Computer Science
Replies
14
Views
4K
Executing Specific Python Lines in VS Code on Ubuntu
    • Programming and Computer Science
Replies
2
Views
855
Are Python strings truly immutable?
    • Programming and Computer Science
Replies
10
Views
2K
Fortran Errors in Fortran code for solving Laplace's equation in 3D
    • Programming and Computer Science
Replies
4
Views
618
Python How Can I Combine Slicing and Reversing to Determine Palindromes in Python?
    • Programming and Computer Science
2
Replies
55
Views
4K
Page not found 404 django-python
    • Programming and Computer Science
Replies
4
Views
3K
Second programming language to get under the belt: python vs C
    • Programming and Computer Science
Replies
9
Views
1K
Re: Python’s Sympy Module and the Cayley-Hamilton Theorem
    • Programming and Computer Science
Replies
3
Views
318
Fortran Fortran Code problems linking to LAPACK library
    • Programming and Computer Science
Replies
2
Views
2K
Python Who can find the largest prime number with the own programmed code?
    • Programming and Computer Science
Replies
22
Views
767

Hot Threads

Recent Insights

Back
Top