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

Click For Summary
The discussion centers around the challenges of translating Fortran code into Python, particularly for a system of 66 ordinary differential equations (ODEs) related to atmospheric emissions. The original poster seeks assistance in running the Fortran code to obtain variable values for comparison with their Python implementation. They express difficulty in understanding the Fortran code, noting that it lacks a main program to call the necessary subroutines and print results. Participants emphasize the need for a main program to execute the subroutines correctly and suggest that without it, the Fortran code cannot be run as intended. There is a consensus that debugging the Fortran code is essential before any translation can be effectively completed. Some participants recommend hiring a Fortran programmer for more specialized help, as the code's current state does not provide the expected output. The conversation highlights the complexities of working with legacy code and the importance of having complete, functional code before attempting to translate it into another programming language.
  • #31
Vick said:

expect​


verb (used with object)
1/ to look forward to; regard as likely to happen; anticipate the occurrence or the coming of
2/ to look for with reason or justification

There are different meanings; my meaning in asking this question is to anticipate the occurrence; whereas you took the second meaning!!
It looks like large sections of your Python code mimic the subroutines in the FORTRAN code. It would help if your Python indicated what it matches with a comment indicating the exact FORTRAN subroutine name. Then someone (not me) might call the FORTRAN subroutines in the same order as your Python code and insert prints that match yours.
 
Technology news on Phys.org
  • #32
Vick said:
Can you please make one and call up these subroutines to print the values of variables RC, DJ, Y, T, DY for TIME=331250?
I don't think anyone here can do that. The Fortran source code file, emep.f, is not a program -- it is essentially a library file that contains 9 or 10 subroutines and a function or two. A main program that exercises these subroutines would need to call them in a particular order and would also need to provide a list of suitably initialized parameters in each call.

The code was revised in October of 2006, but it looks to be plain vanilla Fortran 77 that was cobbled together to do something or other.

In one of the subroutines there is this monstrosity. I've included probably less than 20% of that block of code.

Fortran: 
      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(5
     &7))+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)+R
     &C(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.19
     &D1*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(2
     &0)*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(1
     &60)*Y(44)+0.5D-1*RC(150)*Y(41))+Y(11)*(DJ(6)+DJ(7)+RC(69)*Y(39))+D
     &J(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(1
     &26)*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)
... and similar up through DY(66).

Many of the other subroutines are only slightly more understandable.
 
  • #33
Vick said:
You told me to go to the source!
I told you that it might help to go back to the place where the FORTRAN code you posted originally came from. You already know where that is since you posted a link to the problem page. The people there would be much better able to answer questions of the kind you are asking.
 
  • #34
Vick said:
As the formulation is lengthy, they just pointed out to this Fortran code. source
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.
 
  • #35
FactChecker said:
It would help if your Python indicated what it matches with a comment indicating the exact FORTRAN subroutine name.
Python: 
import numpy as np

'''
C     EMEP MSC-W OZONE MODEL CHEMISTRY
C
C=======================================================================
C
      INTEGER I
C Parameters
      INTEGER NSPEC
      PARAMETER (NSPEC=66)
C     NSPEC : 66    number of species

      FULLNM = 'EMEP problem'
      PROBLM = 'emep'
      TYPE   = 'ODE'
      NEQN   = NSPEC
      NDISC  = 8
      T(0)   =  4D0*3600D0
      T(1)   = 20D0*3600D0
      DO 10 I=1,4
         T(2*I)   = T(0) + 24D0*3600D0*DBLE(I)
         T(2*I+1) = T(1) + 24D0*3600D0*DBLE(I)
   10 CONTINUE
      NUMJAC = .FALSE.
      MLJAC  = NEQN
      MUJAC  = NEQN

      RETURN
      END
'''

# from Fortran code above
#-------------------------
tt = np.empty(10)

tt[0] = 4*3600
tt[1] = 20*3600
for i in range(1,5):
    tt[2*i] = tt[0] +24*3600*i
    tt[2*i+1] = tt[1] +24*3600*i

print(tt)
#-------------------------

'''
      subroutine init(neqn,t,y,yprime,consis)
      integer neqn
      double precision t,y(neqn),yprime(neqn)
      logical consis

      INTEGER I
C
C
C Establishment of initial conditions:
C    completely arbitrary at this stage !
      DO 10 I = 1, 13
          Y(I)=1.D7
 10   CONTINUE
C
      DO 20 I = 14, NEQN
         Y(I)=100.D0
 20   CONTINUE
      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-3
C
      RETURN
      END
'''

# from Fortran code above
#-------------------------
# initial conditions
u = np.empty(66)

for i in range(0,13):
    u[i] = 1.e+7
for i in range(13, 66):
    u[i] = 100.
u[0] = 1.e+9
u[1] = 5.e+9
u[2] = 5.e+9
u[3] = 3.8e+12
u[4] = 3.5e+13
u[13] = 5.e+11
u[37] = 1.e-3
print()
print(u)

#-------------------------

'''
C Parameters
      INTEGER NSPEC, NRC, NDJ
      DOUBLE PRECISION HMIX
      PARAMETER (NSPEC=66, NRC=266, NDJ=16)
      PARAMETER (HMIX=1.2D5)
C     NSPEC : 66    number of species
C     NRC   : 266   size of rate constant array
C     NDJ   : 16    number of photolysis reactions
C     HMIX  : mixing height in cm
C     EMIS1..EMIS13 : emitted species
C
C Listing of species
C     Y(1:NSPEC) =
C        NO,     NO2,    SO2,    CO,     CH4,     C2H6,
C        NC4H10, C2H4,   C3H6,   OXYL,   HCHO,    CH3CHO,
C        MEK,    O3,     HO2,    HNO3,   H2O2,    H2,
C        CH3O2,  C2H5OH, SA,     CH3O2H, C2H5O2,  CH3COO,
C        PAN,    SECC4H, MEKO2,  R2OOH,  ETRO2,   MGLYOX,
C        PRRO2,  GLYOX,  OXYO2,  MAL,    MALO2,   OP,
C        OH,     OD,     NO3,    N2O5,   ISOPRE,  NITRAT,
C        ISRO2,  MVK,    MVKO2,  CH3OH,  RCO3H,   OXYO2H,
C        BURO2H, ETRO2H, PRRO2H, MEKO2H, MALO2H,  MACR,
C        ISNI,   ISRO2H, MARO2,  MAPAN,  CH2CCH3, ISONO3,
C        ISNIR,  MVKO2H, CH2CHR, ISNO3H, ISNIRH,  MARO2H
C
C=======================================================================
C Emissions in molec/(cm**2*s), of NO, NO2, SO2, CO, CH4, C2H6, NC4H10,
C    C2H4, C3H6, O-XYLENE, C2H5OH and ISOPRENE
      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
C
C   distribution of VOC emissions among species:
      PARAMETER (FRAC6=0.07689D0, FRAC7=0.41444D0, FRAC8=0.03642D0,
     +           FRAC9=0.03827D0, FRAC10=0.24537D0, FRAC13=0.13957D0)
      PARAMETER (VEKT6=30.D0, VEKT7=58.D0, VEKT8=28.D0,
     +           VEKT9=42.D0, VEKT10=106.D0, VEKT13=46.D0)
      PARAMETER (VMHC=1.D0/(FRAC6/VEKT6+FRAC7/VEKT7+FRAC8/VEKT8+
     +           FRAC9/VEKT9+FRAC10/VEKT10+FRAC13/VEKT13))
C
C   choose values for NOX and HC emissions in molecules/cm**2xs
      PARAMETER (EMNOX=2.5D11, EMHC=2.5D11)
C
C   rural case
      PARAMETER ( FACISO=1.0D0, FACHC=1.0D0, FACNOX=1.0D0)
C
      PARAMETER (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)
'''# from Fortran code above
#-------------------------
# listing of species

ls = ['NO',     'NO2',    'SO2',    'CO',     'CH4',     'C2H6',
     'NC4H10', 'C2H4',   'C3H6',   'OXYL',   'HCHO',    'CH3CHO',
     'MEK',    'O3',     'HO2',    'HNO3',   'H2O2',    'H2',
     'CH3O2',  'C2H5OH', 'SA',     'CH3O2H', 'C2H5O2',  'CH3COO',
     'PAN',    'SECC4H', 'MEKO2',  'R2OOH',  'ETRO2',   'MGLYOX',
     'PRRO2',  'GLYOX',  'OXYO2',  'MAL',    'MALO2',   'OP',
     'OH',     'OD',     'NO3',    'N2O5',   'ISOPRE',  'NITRAT',
     'ISRO2',  'MVK',    'MVKO2',  'CH3OH',  'RCO3H',   'OXYO2H',
     'BURO2H', 'ETRO2H', 'PRRO2H', 'MEKO2H', 'MALO2H',  'MACR',
     'ISNI',   'ISRO2H', 'MARO2',  'MAPAN',  'CH2CCH3', 'ISONO3',
     'ISNIR',  'MVKO2H', 'CH2CHR', 'ISNO3H', 'ISNIRH',  'MARO2H']# mixing height in cm
HMIX = 1.2e+5

# distribution of VOC emissions among species
FRAC6=0.07689
FRAC7=0.41444
FRAC8=0.03642
FRAC9=0.03827
FRAC10=0.24537
FRAC13=0.13957

VEKT6=30.0
VEKT7=58.0
VEKT8=28.0
VEKT9=42.0
VEKT10=106.0
VEKT13=46.0

VMHC=1.0/(FRAC6/VEKT6+FRAC7/VEKT7+FRAC8/VEKT8+ \
          FRAC9/VEKT9+FRAC10/VEKT10+FRAC13/VEKT13)# values for NOX and HC emissions in molecules/cm**2xs
EMNOX = 2.5e+11
EMHC = 2.5e+11# rural case
FACISO = 1.0
FACHC = 1.0
FACNOX = 1.0

# Emissions in molec/(cm^2*s) of NO, NO2, SO2, CO, CH4, C2H6, NC4H10,
#                                C2H4, C3H6, O-XYLENE, C2H5OH and ISOPRENE
EMIS1 = EMNOX * FACNOX
EMIS2 = 0.0
EMIS3 = EMNOX * FACNOX
EMIS4 = EMHC*10.0 * FACHC
EMIS5 = 0.0
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.5 * FACISO * EMHC
EMIS12= 0.0
EMIS13= EMHC*FRAC13/VEKT13*VMHC* FACHC
#-------------------------

'''
C Compute time-dependent EMEP coefficients
      CALL EMEPCF (TIME, RC, DJ, H2O)
C
C
      M = 2.55D19
      O2 = 5.2D18
      XN2= 1.99D19
C..  pathways for decay of secc4h9o:
C..  Newer assumption, from Atkisnon , 1991
      RPATH3 = 0.65D0
      RPATH4 = 0.35D0
'''# from Fortran code above
#-------------------------
# time-dependent EMEP coefficients
M = 2.55e+19
O2 = 5.2e+18
XN2 = 1.99e+19
RPATH3 = 0.65
RPATH4 = 0.35
#-------------------------
'''
    SUBROUTINE EMEPCF (TIME, RC, DJ, H2O)
      INTEGER NSPEC, NRC, NDJ
      DOUBLE PRECISION TIME, HMIX
      PARAMETER (NSPEC=66, NRC=266, NDJ=16)
      PARAMETER (HMIX=1.2D5)
      DOUBLE PRECISION RC(NRC), DJ(NDJ), H2O
C
C Compute time-dependent EMEP coefficients
C   RC: reaction coefficients
C   DJ: dissociation rate coefficient
C   H2O water vapour concentrations
C
C A and B: DJ=A*exp(-B*SEC)
C    SEC = 1/cos(THETA) where THETA is solar zenith angle
C T temperature in K
C
      DOUBLE PRECISION A(NDJ), B(NDJ), SEC, T

      INTEGER ITIMEH, I24HRS, I
      DOUBLE PRECISION TIMEH, TIMEOD, PI, XLHA, FI, DEKL, XQ, RH, XZ
      DOUBLE PRECISION M, O2, XN2, DELTA

      DATA A/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/
      DATA 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/

      TIMEH=TIME/3600.D0
      ITIMEH=int(TIMEH)
      I24HRS=ITIMEH/24+1
      TIMEOD=TIMEH-(I24HRS-1)*24.D0
C
C Meteorology
C
      PI=4.D0*ATAN(1.0D0)
C   XLHA local hour angle
      XLHA=(1.D0+TIMEOD*3600.D0/4.32D4)*PI
C   FI (Norwegian PHI!) latitude, dekl solar declination
C   here latitude = 50 deg. N
      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))
C   def of temperature variation
C     XP=8.7D-5*TIMEOD*3600.D0-2.83
C     T=8.3D0*SIN(XP)+289.86
C   for simplicity
      T=298.D0
C
C   def of water vapor concentration
      XQ=-7.93D-5*TIMEOD*3600.D0+2.43D0
      RH=23.D0*SIN(XQ)+66.5D0
C
      XZ=(597.3D0-0.57D0*(T-273.16D0))*18.D0/1.986D0*
     1 (1.D0/T-1.D0/273.16D0)
      H2O=6.1078D0*EXP(-XZ)*10.D0/(1.38D-16*T)*RH
'''

# from Fortran code above
#-------------------------
# dissociation rate coefficient
TIME = 331250 # to input
TIMEH=TIME/3600.
ITIMEH=int(TIMEH)
I24HRS=ITIMEH/24+1
TIMEOD=int(24. -(I24HRS-1))+(TIMEH-ITIMEH)
print()
print(TIMEOD)

# XLHA local hour angle
XLHA=(1.0+TIMEOD*3600.0/4.32e+4)*np.pi

# FI (Norwegian PHI!) latitude, dekl solar declination
# here latitude = 50 deg. N
FI=50.0*np.pi/180.0
DEKL=23.5*np.pi/180.0
SEC=1.0/(np.cos(XLHA)*np.cos(FI)*np.cos(DEKL)+np.sin(FI)*np.sin(DEKL))#  def of temperature variation
#     XP=8.7D-5*TIMEOD*3600.D0-2.83
#     T=8.3D0*SIN(XP)+289.86
#   for simplicity
# T temperature in K
T=298.0

#   def of water vapor concentration
XQ=-7.93e-5*TIMEOD*3600.0+2.43
RH=23.0*np.sin(XQ)+66.5

XZ=(597.3-0.57*(T-273.16))*18.0/1.986*(1.0/T-1.0/273.16)
H2O=6.1078*np.exp(-XZ)*10.0/(1.38e-16*T)*RH

A = [1.23e-3,2.00e-4,1.45e-2,2.20e-5,3.00e-6, \
     5.40e-5,6.65e-5,1.35e-5,2.43e-5,5.40e-4, \
     2.16e-4,5.40e-5,3.53e-2,8.94e-2,3.32e-5,2.27e-5]

B = [0.60,   1.40,   0.40,   0.75,   1.25, \
     0.79,   0.60,   0.94,   0.88,   0.79, \
     0.79,   0.79,  0.081,  0.059,   0.57,   0.62]
#-------------------------# SEC = 1./np.cos(THETA) # THETA is solar zenith angle'''
C Calculate  values of photolysis rates DJ(1..16), based
C   upon RGD A & B coefficients and correction factors from HOUGH (1988)
C
      IF (TIMEOD .LT. 4.0D0 .OR. TIMEOD .GE. 20.D0) THEN
C   in the dark:
         DO 100 I = 1,NDJ
            DJ(I)=1.D-40
 100     CONTINUE
      ELSE
C   daytime:
         DO 110 I = 1,NDJ
            DJ(I)=A(I) * EXP( -B(I) * SEC )
            IF( DJ(I) .LT. 0.0D0 ) STOP 'DJ'
 110     CONTINUE
      ENDIF
'''

# from Fortran code above
#-------------------------
# Calculate  values of photolysis rates DJ(1..16), based
#   upon RGD A & B coefficients and correction factors from HOUGH (1988)
DJ = np.empty(16)

if TIMEOD < 4.0 or TIMEOD >= 20.0:
    # in the dark
    for i in range(16):
        DJ[i] = 1.0e-40
else:
    # daytime
    for i in range(16):
        DJ[i] = A[i]*np.exp(-B[i]*SEC)
        if DJ[i] < 0.0:
            break

print()
print(DJ)
#-------------------------

'''
C Set up chemical reaction rate coefficients:
C     16/6/92: inclusion of M, N2, O2 values in rate-constants
C     reaction rate coefficient definition. units: 2-body reactions
C     cm**3/(molecule x s), unimolecular 1.D0/s, 3-body
C     cm**6/(molecule**2 x s)
C
      M  = 2.55D19
      O2 = 5.2D18
      XN2= 1.99D19
C
      DO 120 I = 1,NRC
         RC(I)=0.D0
 120  CONTINUE
c
c
c..A92, assuming 80% N2, 20% O2
      rc(1)  =5.7d-34*(t/300.0D0)**(-2.8D0) * m * o2
c..unchanged, as source of O+NO+O2 reaction rate unknown.
      rc(5)  =9.6d-32*(t/300.0D0)**(-1.6D0) * m
c..estimate from Demore(1990) (=A92) O(d)+N2 and DeMore O(D)+O2
      rc(7)  =2.0d-11*exp(100.D0/t) * m
c.. A92    >>>>>
      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)

.
.
.
to end of Fortran code
'''# from Fortran code above
#-------------------------
# Set up chemical reaction rate coefficients:

RC = np.empty(266)

for i in range(266):
    RC[i] = 0.0

# assuming 80% N2, 20% O2
RC[0]=5.7e-34*(T/300.0)**(-2.8) * M * O2

# unchanged, as source of O+NO+O2 reaction rate unknown.
RC[4]=9.6e-32*(T/300.0)**(-1.6) * M

# estimate from Demore(1990) (=A92) O(d)+N2 and DeMore O(D)+O2
RC[6]  =2.0e-11*np.exp(100.0/T) * M

RC[10] =1.8e-12*np.exp(-1370.0/T)
RC[11] =1.2e-13*np.exp(-2450.0/T)
RC[12] =1.9e-12*np.exp(-1000.0/T)
RC[13] =1.4e-14*np.exp(-600.0/T)
RC[14] =1.8e-11*np.exp(+110.0/T)
RC[16] =3.7e-12*np.exp(240.0/T)

# M.J (from Wayne et al.)
RC[18] =7.2e-14*np.exp(-1414.0/T)
# RC[18] =7.2e-13*np.exp(-1414.0/T)

# M.J suggests that rc(27) be omitted:
#     RC[27] =8.5e-13*np.exp(-2450.0/T)
#  My change to get similar to A92 troe.
RC[28] =7.1e+14*np.exp(-11080.0/T)

RC[29] =4.8e-11*np.exp(+250.0/T)

# De More,1990 .. no change : oh + h2o2
RC[30] =2.9e-12*np.exp(-160.0/T)

# : oh + h2
RC[32] =7.7e-12*np.exp(-2100.0/T)

# My, similar to DeMore et al complex : oh+hno3
RC[34] =1.0e-14*np.exp(785.0/T)

# Mike Jenkin`s suggestion for ho2 + ho2 reactions: (from DeMore et al.)
RC[35] =2.3e-13*np.exp(600.0/T)
RC[35] = RC[36] + M * 1.7e-33*np.exp(1000.0/T)
RC[35] = RC[35] * (1.0 + 1.4e-21 * H2O *np.exp(2200.0/T))

RC[58] =3.9e-12*np.exp(-1885.0/T)
# : ch3o2 + no
RC[59] =4.2e-12*np.exp(180.0/T)
# A92 + A90 assumption that ka = 0.5D0 * k
RC[60] =5.5e-14*np.exp(365.0/T)

# A92 + A90 assumption that kb = 0.5D0 * k
RC[61] =5.5e-14*np.exp(365.0/T)

RC[62] =3.3e-12*np.exp(-380.0/T)
# : ho2 + ch3o2
RC[64] =3.8e-13*np.exp(780.0/T)

# new: ch3ooh + oh -> hcho + oh
RC[66] =1.0e-12*np.exp(190.0/T)
# new: ch3ooh + oh -> ch3o2
RC[67] =1.9e-12*np.exp(190.0/T)

RC[70] =7.8e-12*np.exp(-1020.0/T)
# new: c2h5o2 + ho2 -> c2h5ooh (r2ooh)
RC[73] =6.5e-13*np.exp(650.0/T)

RC[74] =5.6e-12*np.exp(310.0/T)
# TOR90 assumption w.r.t. rc(67) : c2h5ooh + oh -> ch3cho + oh
#     RC[75] = 5.8 * RC[66]
RC[75] = 5.8e-12*np.exp(190.0/T)
# : approximation to troe expression
RC[77] =1.34e+16*np.exp(-13330.0/T)
# additional reactions :-
# : ho2 + ch3coo2 -> rco3h
RC[87] =1.3e-13*np.exp(1040.0/T)
# : ho2 + ch3coo2 -> rco2h + o3
RC[88] =3.0e-13*np.exp(1040.0/T)

RC[93] =2.8e-12*np.exp(530.0/T)
# D & J, gives results very close to A92
RC[80] =1.64e-11*np.exp(-559.0/T)RC[82] =RC[59]
RC[104]=RC[59]

RC[109]=RC[59]
# A90/PS  isnir + no
RC[161]=RC[59]
# A90/PS  isono3 + no
RC[163]=RC[59]

# From RGD, but very similar to A92 Troe
RC[108]=1.66e-12*np.exp(474.0/T)

RC[111]=1.2e-14*np.exp(-2630.0/T)

RC[122]=6.5e-15*np.exp(-1880.0/T)

RC[125] = RC[59]

RC[219] = RC[59]

RC[235] = RC[59]
#  natural voc reactions #
#  isoprene + o3 -> products
RC[149] = 12.3e-15*np.exp(-2013.0/T)
#  isoprene + oh -> isro2
RC[150] = 2.54e-11*np.exp(410.0/T)
#  isoprene-RO2 + no
RC[151] = RC[59]
#  methylvinylketone (mvk) + oh
RC[152] = 4.13e-12*np.exp(452.0/T)
#  mvko2 + no
RC[153] = RC[59]
#  macr + oh
RC[157] = 1.86e-11*np.exp(175.0/T)
#  ch2cch3 + no
RC[158] = RC[59]
# mvk + o3
RC[159] = 4.32e-15*np.exp(-2016.0/T)

#     RC[254] = 9.9e-16*np.exp(-731.0/T)

#     aerosol formation and depositions ..x
#     parameterization of heterogeneous loss in 1./s for humidities
#     less than 90%.  applied to hno3, h2o2, h2so4, ch3o2h.
RC[42]=5.e-6
if(RH > 0.90):
    RC[42]=1.e-4

RC[43] =RC[42]
RC[44] =RC[42]RC[7]  =2.2e-10
# My (?) , similar to A92 troe expression.
RC[19] =1.4e-12
# My (?) , similar to A92 troe expression.
RC[20] =1.4e-11
# RGD Note (=DeMore?)
RC[25] =4.1e-16
# RGD
RC[38] =1.35e-12
# RGD
RC[39] =4.0e-17
# A92, assuming all products are ch3cho
RC[63] =3.2e-12
# A92, with temperature dependance neglected.
RC[65] =9.6e-12

RC[68] =5.8e-16
RC[69] =2.4e-13
RC[71] =8.9e-12

# A92 : approximation to troe expression
RC[76] =1.0e-11
# : ch3coo2 + no
RC[78] =2.0e-11
#: sum of two pathways
RC[79] =1.1e-11
# cmj   RC[83] =2.5e-14
#   kho2ro2 : estimate of rate of ho2 + higher RO2, suggested by Yvonne
RC[84] =1.0e-11
# ignoring slight temperature dependance
RC[85] =1.15e-12
# MJ suggestion.. combine oh + secc4h9o2, meko2 rates   =>
#     RC[86]= RC[67] + RC[85], approx. at 298
RC[86]= 4.8e-12
# new A92 : oh + ch3co2h -> ch3o2
RC[89] =8.0e-13
# Approximates to A92 Troe ...
RC[124]=2.86e-11
# rate for ch2chr + oh, = k68+k125 (propene), Yv.
RC[145] = 3.2e-11
# rate for isono3h + oh, = k156 (isro2h), Yv.
RC[146] = 2.0e-11
# MY GUESS rate of oh + mvko2h
RC[147] = 2.2e-11
# MY GUESS rate of oh + other biogenic ooh
RC[148] = 3.7e-11
#   kho2ro2 : estimate of rate of ho2 + higher RO2, suggested by Yvonne
#  isro2 + ho2
RC[154] =RC[84]
#  isro2h + oh
RC[155] =2.0e-11
#   isro2h + o3
RC[156] =8.0e-18
#   isni + oh
RC[160] =3.35e-11
#  isopre + no3
RC[162] =7.8e-13
#  RC[162] =7.8e-16
# Unchanged, also in IVL scheme
RC[218]=2.0e-11

RC[220]=1.1e-11

RC[221]=1.7e-11
# MJ suggestion.. combine oh + malo2h rates   =>
#     RC[222]= RC[67] + RC[218]
RC[222]= 2.4e-11

RC[233]=1.37e-11
# MJ suggestion.. combine rc(68) with rc(234) =>
#     RC[234]= RC[67] + RC[233]
RC[234]= 1.7e-11

#         deposition loss rate coefficients vd/hmix, vd in cm/s.
#         hno3     calculated     rc(46)
#         so2      0.5            rc(47)
#         h2o2     0.5            rc(47)
#         no2      0.2            rc(48)
#         o3       0.5            rc(49)
#         pan      0.2            rc(50)
#         h2so4    0.1            rc(51)
# simple approx. for now - reduce all vg by 4 at night to allow
#    for surface inversion......
delta=1.0
if(TIMEOD >= 20.0 or TIMEOD < 4.0):
    delta=0.25
#         if(timeod.ge.20.D0.or.timeod.le.4.D0) delta=0.25D0
#         if (night) delta=0.25D0
RC[45] =2.0 * delta /HMIX
RC[46] =0.5 * delta /HMIX
RC[47] =0.2 * delta /HMIX
RC[48] =0.5 * delta /HMIX
RC[49] =0.2 * delta /HMIX
RC[50] =0.1 * delta /HMIX
# dep. of organic peroxides = 0.5 cms-1
RC[51] = 0.5 *delta /HMIX
# dep. of ketones, RCHO  = 0.3 cms-1
RC[52] = 0.3 *delta /HMIX

print()
print(RC)
#-------------------------'''
C=======================================================================
C
C
      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(5
     &7))+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)+R
     &C(79)*(Y(24)+Y(57))+0.2D1*RC(15)*Y(39))+RC(60)*Y(1)*(Y(19)+Y(26)+Y

etc

'''

# from Fortran code above
#-------------------------
# ODEs
Du = np.empty(66)
Du[1-1] = DJ[3-1]*u[2-1]+DJ[13-1]*u[39-1]+RC[19-1]*u[2-1]*u[39-1]+EMIS1/HMIX- \
        (RC[5-1]*u[36-1]+RC[11-1]*u[14-1]+RC[17-1]*u[15-1]+RC[72-1]*u[23-1]+RC[79-1]* \
         (u[24-1]+ u[57-1])+RC[15-1]*u[39-1]+RC[60-1]*(u[19-1]+u[26-1]+u[27-1]+u[29-1]+u[31-1]+u[33-1]+ \
                                             u[35-1]+u[43-1]+u[45-1]+u[59-1]+u[61-1]+u[60-1]))*u[1-1]

Du[2-1] = u[1-1]*(RC[5-1]*u[36-1]+RC[11-1]*u[14-1]+RC[17-1]*u[15-1]+RC[72-1]*u[23-1]+ \
              RC[79-1]*(u[24-1]+u[57-1])+0.2e+1*RC[15-1]*u[39-1])+RC[60-1]*u[1-1]* \
              (u[19-1]+u[26-1]+u[27-1]+u[29-1]+u[31-1]+u[33-1]+u[35-1]+u[59-1])+ \
              RC[60-1]*u[1-1]*(0.86*u[43-1]+0.19e+1*u[61-1]+0.11e+1*u[60-1]+0.95*u[45-1])+ \
              DJ[14-1]*u[39-1]+DJ[5-1]*u[16-1]+DJ[15-1]*u[40-1]+RC[29-1]*u[40-1]+RC[78-1]* \
              (u[25-1]+u[58-1])-(DJ[3-1]+RC[12-1]*u[14-1]+RC[20-1]*u[39-1]+RC[21-1]*u[37-1]+RC[48-1]+RC[77-1]* \
                             (u[24-1]+u[57-1]))*u[2-1]

Du[3-1] = EMIS3/HMIX-(RC[39-1]*u[37-1]+RC[40-1]*u[19-1]+RC[47-1])*u[3-1]

Du[4-1] = EMIS4/HMIX+u[37-1]*(RC[66-1]*u[11-1]+2.0*RC[221-1]*u[32-1]+RC[222-1]* u[30-1])+ \
        u[14-1]*(0.44*RC[112-1]*u[8-1]+0.4*RC[123-1]*u[9-1]+0.5e-1*RC[160-1]*u[44-1]+0.5e-1* \
               RC[150-1]*u[41-1])+u[11-1]*(DJ[6-1]+DJ[7-1]+RC[69-1]*u[39-1])+DJ[8-1]*u[12-1]+DJ[11-1]*u[30-1]+2.0*DJ[7-1]*u[32-1]-RC[70-1]*u[37-1]*u[4-1]

Du[5-1] = EMIS5/HMIX+0.7e-1*RC[123-1]*u[14-1]*u[9-1]-RC[59-1]*u[37-1]*u[5-1]

Du[6-1] = EMIS6/HMIX-RC[71-1]*u[37-1]*u[6-1]

Du[7-1] = EMIS7/HMIX-RC[81-1]*u[37-1]*u[7-1]

Du[8-1] = EMIS8/HMIX-(RC[109-1]*u[37-1]+RC[112-1]*u[14-1])*u[8-1]

Du[9-1] = EMIS9/HMIX+0.7e-1*RC[150-1]*u[14-1]*u[41-1]-(RC[123-1]*u[14-1]+RC[125-1]*u[37-1])*u[9-1]

Du[10-1] = EMIS10/HMIX-RC[234-1]*u[37-1]*u[10-1]

Du[11-1] = u[19-1]*(RC[60-1]*u[1-1]+(2.0*RC[61-1]+RC[62-1])*u[19-1]+RC[80-1]*u[24-1]+RC[40-1]*u[3-1])+ \
         u[37-1]*(RC[63-1]*u[46-1]+RC[67-1]*u[22-1])+u[1-1]*RC[60-1]* \
         (2.0*u[29-1]+u[31-1]+0.74*u[43-1]+0.266*u[45-1]+0.15*u[60-1])+u[14-1]* \
         (0.5*RC[123-1]*u[9-1]+RC[112-1]*u[8-1]+0.7*RC[157-1]*u[56-1]+0.8*RC[160-1]*u[44-1]+0.8*RC[150-1]* \
          u[41-1])+2.0*DJ[7-1]*u[32-1]+DJ[16-1]* \
          (u[22-1]+0.156e+1*u[50-1]+u[51-1])-(RC[66-1]*u[37-1]+DJ[6-1]+DJ[7-1]+RC[69-1]*u[39-1]+RC[53-1])*u[11-1]

Du[12-1] = u[1-1]*(RC[72-1]*u[23-1]+RC[83-1]*u[26-1]*RPATH4+RC[105-1]*u[27-1]+RC[126-1]*u[31-1]+ \
               0.95*RC[162-1]*u[61-1]+0.684*RC[154-1]*u[45-1])+u[37-1]* \
               (RC[64-1]*u[20-1]+RC[76-1]*u[28-1]+RC[76-1]*u[50-1])+0.5*RC[123-1]*u[14-1]*u[9-1]+0.4e-1* \
               RC[160-1]*u[14-1]*u[44-1]+DJ[16-1]*(u[28-1]+0.22*u[50-1]+0.35*u[49-1]+u[51-1]+u[52-1])- \
               (DJ[8-1]+RC[75-1]*u[37-1]+RC[53-1])*u[12-1]

Du[13-1] = RC[83-1]*u[1-1]*u[26-1]*RPATH3+(0.65*DJ[16-1]+RC[76-1]*u[37-1])*u[49-1]+RC[76-1]*u[37-1]*u[51-1]+ \
         (RC[159-1]*u[1-1]+DJ[16-1])*u[59-1]+0.95*RC[162-1]*u[1-1]*u[61-1]-(DJ[9-1]+RC[86-1]*u[37-1]+RC[53-1])*u[13-1]

Du[14-1] = RC[1-1]*u[36-1]+RC[89-1]*u[15-1]*u[24-1]-(RC[11-1]*u[1-1]+RC[12-1]*u[2-1]+RC[13-1]*u[37-1]+RC[14-1]*u[15-1]+RC[49-1]+ \
                                         RC[112-1]*u[8-1]+RC[123-1]*u[9-1]+RC[157-1]*u[56-1]+RC[160-1]*u[44-1]+RC[150-1]* \
                                         u[41-1]+DJ[1-1]+DJ[2-1])*u[14-1]

s1 = u[37-1]*(RC[13-1]*u[14-1]+RC[31-1]*u[17-1]+RC[33-1]*u[18-1]+RC[39-1]*u[3-1]+RC[63-1]*u[46-1]+RC[64-1]* \
            u[20-1]+RC[66-1]*u[11-1]+RC[70-1]*u[4-1]+RC[221-1]*u[32-1])+u[19-1]* \
            (RC[40-1]*u[3-1]+2.0*RC[61-1]*u[19-1]+0.5*RC[80-1]*u[24-1])+DJ[11-1]*u[30-1]+u[1-1]*RC[60-1]* \
            (u[19-1]+u[29-1]+u[31-1]+u[33-1]+u[35-1]+0.95*u[45-1]+u[26-1]*RPATH3+0.78*u[43-1]+u[59-1]+ \
             0.5e-1*u[61-1]+0.8*u[60-1])+RC[72-1]*u[1-1]*u[23-1]+DJ[8-1]*u[12-1]

Du[15-1] = s1+2.0*DJ[6-1]*u[11-1]+DJ[16-1]*(u[22-1]+u[28-1]+0.65*u[49-1]+u[50-1]+u[51-1]+u[48-1]+u[53-1])+ \
         u[39-1]*(RC[26-1]*u[17-1]+RC[69-1]*u[11-1])+u[14-1]*(0.12*RC[112-1]*u[8-1]+0.28*RC[123-1]*u[9-1]+ \
                                                  0.6e-1*RC[160-1]*u[44-1])+0.6e-1*RC[150-1]*u[14-1]*u[41-1]- \
                                                  (RC[14-1]*u[14-1]+RC[17-1]*u[1-1]+RC[30-1]*u[37-1]+2.0*RC[36-1]* \
                                                   u[15-1]+RC[65-1]*u[19-1]+RC[74-1]*u[23-1]+(RC[88-1]+RC[89-1])*u[24-1]+ \
                                                   RC[85-1]*(u[26-1]+u[29-1]+u[31-1]+u[27-1]+u[57-1]+u[45-1]+u[61-1]+u[59-1]+u[33-1]+ \
                                                           u[35-1]+u[43-1]+u[60-1]))*u[15-1]

Du[16-1] = RC[21-1]*u[2-1]*u[37-1]+u[39-1]*(RC[26-1]*u[17-1]+RC[69-1]*u[11-1])-(RC[35-1]*u[37-1]+DJ[5-1]+RC[45-1])*u[16-1]

Du[17-1] = RC[36-1]*u[15-1]**2-(RC[31-1]*u[37-1]+DJ[4-1]+RC[43-1]+RC[26-1]*u[39-1]+RC[47-1])*u[17-1]

Du[18-1] = DJ[7-1]*u[11-1]+u[14-1]*(0.13*RC[112-1]*u[8-1]+0.7e-1*RC[123-1]*u[9-1])-RC[33-1]*u[37-1]*u[18-1]

Du[19-1] = u[37-1]*(RC[59-1]*u[5-1]+RC[68-1]*u[22-1])+u[24-1]*(RC[79-1]*u[1-1]+2.0*RC[94-1]*u[24-1])+ \
         DJ[8-1]*u[12-1]+DJ[16-1]*u[47-1]+0.31*RC[123-1]*u[14-1]*u[9-1]-(RC[40-1]*u[3-1]+RC[60-1]*u[1-1]+ \
                                                           2.0*RC[61-1]*u[19-1]+2.0*RC[62-1]*u[19-1]+ \
                                                           RC[65-1]*u[15-1]+0.5*RC[80-1]*u[24-1])*u[19-1]

Du[20-1] = EMIS13/HMIX-RC[64-1]*u[37-1]*u[20-1]

Du[21-1] = (RC[40-1]*u[19-1]+RC[39-1]*u[37-1])*u[3-1]+0.5e-1*EMIS3/HMIX-RC[51-1]*u[21-1]

Du[22-1] = RC[65-1]*u[15-1]*u[19-1]-(RC[43-1]+DJ[16-1]+(RC[67-1]+RC[68-1])*u[37-1])*u[22-1]

Du[23-1] = u[37-1]*(RC[71-1]*u[6-1]+RC[68-1]*u[28-1])+0.35*DJ[16-1]*u[49-1]+RC[83-1]*u[1-1]*u[26-1]* \
         RPATH4+DJ[9-1]*u[13-1]-(RC[72-1]*u[1-1]+RC[74-1]*u[15-1])*u[23-1]

Du[24-1] = u[37-1]*(RC[75-1]*u[12-1]+RC[222-1]*u[30-1]+RC[68-1]*u[47-1])+RC[105-1]* \
         u[1-1]*u[27-1]+RC[78-1]*u[25-1]+DJ[11-1]*u[30-1]+DJ[9-1]*u[13-1]+DJ[16-1]*u[52-1]+ \
         0.684*RC[154-1]*u[1-1]*u[45-1]-(RC[77-1]*u[2-1]+RC[79-1]*u[1-1]+RC[80-1]*u[19-1]+ \
                                   2.0*RC[94-1]*u[24-1]+(RC[88-1]+RC[89-1])*u[15-1])*u[24-1]

Du[25-1] = RC[77-1]*u[24-1]*u[2-1]-(RC[50-1]+RC[78-1])*u[25-1]

Du[26-1] = u[37-1]*(RC[81-1]*u[7-1]+RC[68-1]*u[49-1])-(RC[83-1]*u[1-1]+RC[85-1]*u[15-1])*u[26-1]

Du[27-1] = u[37-1]*(RC[86-1]*u[13-1]+RC[87-1]*u[52-1])-(RC[105-1]*u[1-1]+RC[85-1]*u[15-1])*u[27-1]

Du[28-1] = RC[74-1]*u[15-1]*u[23-1]-((RC[76-1]+RC[68-1])*u[37-1]+DJ[16-1]+RC[52-1])*u[28-1]

Du[29-1] = u[37-1]*(RC[109-1]*u[8-1]+RC[68-1]*u[50-1])-(RC[110-1]*u[1-1]+RC[85-1]*u[15-1])*u[29-1]

Du[30-1] = RC[236-1]*u[1-1]*u[33-1]+RC[220-1]*u[1-1]*u[35-1]+0.266*RC[154-1]*u[1-1]*u[45-1]+ \
         0.82*RC[160-1]*u[14-1]*u[44-1]+DJ[16-1]*(u[48-1]+u[53-1])-(DJ[11-1]+RC[222-1]*u[37-1])*u[30-1]

Du[31-1] = u[37-1]*(RC[125-1]*u[9-1]+RC[68-1]*u[51-1])-(RC[126-1]*u[1-1]+RC[85-1]*u[15-1])*u[31-1]

Du[32-1] = RC[220-1]*u[1-1]*u[35-1]+DJ[16-1]*u[53-1]-(2.0*DJ[7-1]+RC[221-1]*u[37-1])*u[32-1]

Du[33-1] = u[37-1]*(RC[234-1]*u[10-1]+RC[235-1]*u[48-1])-(RC[236-1]*u[1-1]+RC[85-1]*u[15-1])*u[33-1]

Du[34-1] = RC[236-1]*u[1-1]*u[33-1]+DJ[16-1]*u[48-1]-RC[219-1]*u[37-1]*u[34-1]

Du[35-1] = u[37-1]*(RC[219-1]*u[34-1]+RC[223-1]*u[53-1])-(RC[220-1]*u[1-1]+RC[85-1]*u[15-1])*u[35-1]

Du[36-1] = DJ[1-1]*u[14-1]+DJ[3-1]*u[2-1]+DJ[14-1]*u[39-1]+RC[7-1]*u[38-1]+0.2*RC[160-1]*u[14-1]* \
         u[44-1]+0.3*RC[150-1]*u[14-1]*u[41-1]-(RC[1-1]+RC[5-1]*u[1-1])*u[36-1]s1 = 2.0*RC[8-1]*H2O*u[38-1]+u[15-1]*(RC[14-1]*u[14-1]+RC[17-1]*u[1-1])+2.0*DJ[4-1]*u[17-1]+DJ[5-1]*u[16-1]
s2 = s1+DJ[16-1]*(u[22-1]+u[28-1]+u[47-1]+u[49-1]+u[50-1]+u[52-1]+u[48-1]+u[53-1])
s3 = s2+u[14-1]*(0.15*RC[123-1]*u[9-1]+0.8e-1*RC[160-1]*u[44-1])
s4 = s3
s6 = 0.55*RC[150-1]*u[14-1]*u[41-1]
s8 = -1
s11 = RC[222-1]*u[30-1]+RC[75-1]*u[12-1]+RC[81-1]*u[7-1]+RC[87-1]*u[52-1]+RC[86-1]* \
      u[13-1]+RC[235-1]*u[48-1]+RC[109-1]*u[8-1]+RC[125-1]*u[9-1]+RC[234-1]*u[10-1]+ \
      RC[223-1]*u[53-1]+RC[219-1]*u[34-1]+RC[31-1]*u[17-1]+RC[21-1]*u[2-1]+RC[148-1]* \
      u[62-1]+RC[64-1]*u[20-1]+RC[39-1]*u[3-1]+RC[71-1]*u[6-1]
s10 = s11+RC[30-1]*u[15-1]+RC[59-1]*u[5-1]+RC[70-1]*u[4-1]+RC[35-1]*u[16-1]+RC[13-1]* \
      u[14-1]+RC[221-1]*u[32-1]+RC[68-1]*(u[22-1]+u[28-1]+u[47-1]+u[50-1]+u[51-1]+u[49-1])+ \
      RC[66-1]*u[11-1]+RC[151-1]*u[41-1]+RC[153-1]*u[44-1]+RC[63-1]*u[46-1]+RC[33-1]*u[18-1]+ \
      RC[158-1]*u[54-1]+RC[146-1]*u[63-1]+RC[149-1]*(u[65-1]+u[66-1])+RC[147-1]*u[64-1]+ \
      RC[161-1]*u[55-1]
s11 = u[37-1]
s9 = s10*s11
s7 = s8*s9
s5 = s6+s7

Du[37-1] = s4+s5

Du[38-1] = DJ[2-1]*u[14-1]-(RC[7-1]+RC[8-1]*H2O)*u[38-1]

Du[39-1] = (RC[29-1]+DJ[15-1])*u[40-1]+RC[12-1]*u[14-1]*u[2-1]+RC[35-1]*u[37-1]*u[16-1]- \
         (RC[15-1]*u[1-1]+RC[26-1]*u[17-1]+RC[163-1]*u[41-1]+RC[19-1]*u[2-1]+RC[20-1]* \
          u[2-1]+DJ[13-1]+DJ[14-1]+RC[69-1]*u[11-1])*u[39-1]

Du[40-1] = RC[20-1]*u[39-1]*u[2-1]-(RC[29-1]+DJ[15-1]+RC[45-1])*u[40-1]

Du[41-1] = EMIS11/HMIX-(RC[151-1]*u[37-1]+RC[163-1]*u[39-1]+RC[150-1]*u[14-1])*u[41-1]

Du[42-1] = RC[45-1]*u[16-1]+2.0*RC[44-1]*u[40-1]-RC[51-1]*u[42-1]

Du[43-1] = u[37-1]*(RC[151-1]*u[41-1]+RC[156-1]*u[56-1])+0.12*RC[152-1]*u[1-1]* \
         u[43-1]-(RC[152-1]*u[1-1]+RC[155-1]*u[15-1])*u[43-1]

Du[44-1] = RC[60-1]*u[1-1]*(0.42*u[43-1]+0.5e-1*u[60-1])+0.26*RC[150-1]*u[14-1]* \
         u[41-1]-(RC[153-1]*u[37-1]+RC[160-1]*u[14-1])*u[44-1]

Du[45-1] = RC[153-1]*u[44-1]*u[37-1]+RC[148-1]*u[37-1]*u[62-1]-(RC[154-1]*u[1-1]+RC[85-1]*u[15-1])*u[45-1]

Du[46-1] = RC[62-1]*u[19-1]**2-RC[63-1]*u[37-1]*u[46-1]

Du[47-1] = RC[88-1]*u[15-1]*u[24-1]-(RC[68-1]*u[37-1]+DJ[16-1]+RC[52-1])*u[47-1]

Du[48-1] = RC[85-1]*u[15-1]*u[33-1]-(RC[235-1]*u[37-1]+DJ[16-1]+RC[52-1])*u[48-1]

Du[49-1] = RC[85-1]*u[15-1]*u[26-1]-((RC[76-1]+RC[68-1])*u[37-1]+DJ[16-1]+RC[52-1])*u[49-1]

Du[50-1] = RC[85-1]*u[15-1]*u[29-1]-((RC[76-1]+RC[68-1])*u[37-1]+DJ[16-1]+RC[52-1])*u[50-1]

Du[51-1] = RC[85-1]*u[15-1]*u[31-1]-((RC[76-1]+RC[68-1])*u[37-1]+DJ[16-1]+RC[52-1])*u[51-1]

Du[52-1] = RC[85-1]*u[15-1]*u[27-1]-(RC[87-1]*u[37-1]+DJ[16-1]+RC[52-1])*u[52-1]

Du[53-1] = RC[85-1]*u[15-1]*u[35-1]-(RC[223-1]*u[37-1]+DJ[16-1]+RC[52-1])*u[53-1]

Du[54-1] = RC[60-1]*u[1-1]*(0.32*u[43-1]+0.1*u[60-1])+0.67*RC[150-1]*u[14-1]*u[41-1]-RC[158-1]*u[37-1]*u[54-1]

Du[55-1] = RC[60-1]*u[1-1]*(0.14*u[43-1]+0.5e-1*u[45-1]+0.85*u[60-1])-RC[161-1]*u[37-1]*u[55-1]

Du[56-1] = RC[155-1]*u[15-1]*u[43-1]-(RC[156-1]*u[37-1]+RC[157-1]*u[14-1]+RC[52-1])*u[56-1]

Du[57-1] = 0.5*RC[158-1]*u[37-1]*u[54-1]+RC[78-1]*u[58-1]+RC[149-1]*u[37-1]*u[66-1]- \
         (RC[77-1]*u[2-1]+RC[79-1]*u[1-1]+RC[85-1]*u[15-1])*u[57-1]

Du[58-1] = RC[77-1]*u[57-1]*u[2-1]-(RC[50-1]+RC[78-1])*u[58-1]

Du[59-1] = RC[79-1]*u[1-1]*u[57-1]+RC[146-1]*u[37-1]*u[63-1]-(RC[159-1]*u[1-1]+RC[85-1]*u[15-1])*u[59-1]

Du[60-1] = RC[163-1]*u[39-1]*u[41-1]+RC[147-1]*u[37-1]*u[64-1]-(RC[164-1]*u[1-1]+RC[85-1]*u[15-1])*u[60-1]

Du[61-1] = RC[161-1]*u[37-1]*u[55-1]+RC[149-1]*u[37-1]*u[65-1]-(RC[162-1]*u[1-1]+RC[85-1]*u[15-1])*u[61-1]

Du[62-1] = RC[85-1]*u[15-1]*u[45-1]-(RC[148-1]*u[37-1]+RC[52-1])*u[62-1]

Du[63-1] = RC[85-1]*u[15-1]*u[59-1]-(RC[146-1]*u[37-1]+RC[52-1])*u[63-1]

Du[64-1] = RC[85-1]*u[15-1]*u[60-1]-(RC[147-1]*u[37-1]+RC[52-1])*u[64-1]

Du[65-1] = RC[85-1]*u[15-1]*u[61-1]-(RC[149-1]*u[37-1]+RC[52-1])*u[65-1]

Du[66-1] = RC[85-1]*u[15-1]*u[57-1]-(RC[149-1]*u[37-1]+RC[52-1])*u[66-1]

print()
print(Du)
#-------------------------
 
  • #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.
 
  • #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.
 
  • #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:
Reply
  • 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.)
 
Reply
  • Like
Likes member 657093
  • #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.
 
Reply
  • Like
Likes member 657093
  • #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.
 
Reply
  • 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.
 
Reply
  • Like
Likes member 657093
  • #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...
 
  • #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.
 
Reply
  • Like
Likes FactChecker
  • #59
Moderator's note: Thread reopened to allow @Vick to post an update.
 
Reply
  • Like
Likes member 657093
  • #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
 
Reply
  • Like
Likes Filip Larsen, pbuk and FactChecker

Similar threads

Trying to install software on a new PC to run Python code I wrote in 2021
  • · Replies 3 ·
Replies
3
Views
2K
Python Porting Python Code to Fortran for Parallel Computing
  • · Replies 14 ·
Replies
14
Views
5K
[SOLVED] Python - Issues Taking an Input as an Integer (or Converting It)
  • · Replies 1 ·
Replies
1
Views
2K
Python Executing Specific Python Lines in VS Code on Ubuntu
  • · Replies 2 ·
Replies
2
Views
1K
Python Are Python strings truly immutable?
  • · Replies 10 ·
Replies
10
Views
3K
Python How Can I Combine Slicing and Reversing to Determine Palindromes in Python?
  • · Replies 55 ·
2
Replies
55
Views
6K
Python Page not found 404 django-python
  • · Replies 4 ·
Replies
4
Views
5K
Python Second programming language to get under the belt: python vs C
  • · Replies 9 ·
Replies
9
Views
2K
Python Re: Python’s Sympy Module and the Cayley-Hamilton Theorem
  • · Replies 3 ·
Replies
3
Views
1K
Trying to run a big python job with AWS (EC2), is there a better way?
  • · Replies 7 ·
Replies
7
Views
2K