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

[member 657093]

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]

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]

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]

This is turning into r/choosingbeggars

 

[PeterDonis]

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]

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]

Thread closed temporarily for Moderation...

 

[berkeman]

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]

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

 

[member 657093]

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:

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

 

[member 657093]

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

 

[BvU]

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 ?

$\$

 

[member 657093]

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

 

[harborsparrow]

I will just say that in my experience the biggest hazard in translating numerical computations from one language to another may have to do with differences in how the two languages store and organize arrays internally. I don't know what Python does but there is a good chance it is different from what Fortran does.

 

[D H]

I will just say that in my experience the biggest hazard in translating numerical computations from one language to another may have to do with differences in how the two languages store and organize arrays internally. I don't know what Python does but there is a good chance it is different from what Fortran does.

This is a hazard in at least two ways: indexing (is the initial index one or zero), and ordering (row major vs column major). Fortran and MATLAB use one-based indexing and column major order. C, C++, and several other languages use zero-based indexing and row major order. The zero vs one based indexing can be overcome. The order issue is harder to overcome.

As an example, the initial versions of the book "Numerical Recipes in C" were a direct transliteration of the original FORTRAN code (not Fortran, but FORTRAN. FORTRAN was the original name, which means very old code). They used some clever (or perhaps not clever) hacks to get around the indexing issue. They did not do anything with regard to multidimensional arrays, which meant that the memory-efficient FORTRAN algorithms became memory-inefficient in C. Recoding algorithms that use multidimensional arrays is non-trivial. Transliteration oftentimes results in poor performing code.

As far as Python goes, numpy does offer an "order=C" (the default) and "order=F" option on creating and populating a multidimensional array, which might offer some improvements. However, the slowness of Python vs C, C++, Fortran, and Julia typically removes all of those improvements.