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

[member 657093]

TL;DR Summary

Debug Fortran Code

Hi,

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.

 

[Filip Larsen]

Perhaps you can use https://www.onlinegdb.com/online_fortran_compiler or similar.

 

[member 657093]

Perhaps you can use https://www.onlinegdb.com/online_fortran_compiler or similar.

I don't know Fortran even though I've already tried this link. But I can't get it to give me the values for my variables!!

 

[FactChecker]

I don't see anywhere that the final solution is printed. It looks like you need to write a main program that calls these subroutines and then prints the values. You should see if the reference that had this code also has a sample of what and how to call the subroutines and print results.

 

[member 657093]

I don't see anywhere that the final solution is printed. It looks like you need to write a main program that calls these subroutines and then prints the values. You should see if the reference that had this code also has a sample of what and how to call the subroutines and print results.

This piece of code is the actual formulation of a system of 66 ODEs (66 species) to solve for atmospheric emissions of several compounds. I am trying to replicate the formulation in Python because my solver is Python-based. However, if my formulation is wrong then I won't be able to solve it. As I don't understand Fortran, I need someone to provide me with the values of the variables in the actual formulation like RC, DJ, Y, T for TIME=331250 for comparison with my Python formulation.

 

[FactChecker]

This piece of code is the actual formulation of a system of 66 ODEs (66 species) to solve for atmospheric emissions of several compounds. I am trying to replicate the formulation in Python because my solver is Python-based. However, if my formulation is wrong then I won't be able to solve it. As I don't understand Fortran, I need someone to provide me with the values of the variables in the actual formulation like RC, DJ, Y, T for TIME=331250 for comparison with my Python formulation.

Right, but something has to call these subroutines in the correct way to get answers. An example of that would really help. Someone might know FORTRAN very well (Not me. I am rusty and don't have a compiler.) and still not know how to use these subroutines.
If you know how to do it, even with your Python version, that would probably help a lot.

 

[member 657093]

Right, but something has to call these subroutines in the correct way to get answers. An example of that would really help. Someone might know FORTRAN very well (Not me. I am rusty and don't have a compiler.) and still not know how to use these subroutines.
If you know how to do it, even with your Python version, that would probably help a lot.

You want my Python formulation?

 

[PeterDonis]

You want my Python formulation?

Please do not post code in attachments. You can use BBCode tags to put code directly in your post.

 

[member 657093]

Please do not post code in attachments. You can use BBCode tags to put code directly in your post.

What if the code is very long?

 

[PeterDonis]

What if the code is very long?

It will still be more readable in a BBCode code block than an attachment.

Or, if you have it online already in some accessible location (like a public github repository), you can post a link to that.

Also, if the code is really that long and complicated, expecting people to help you rewrite it in a forum discussion here might be too much to expect.

 

[member 657093]

It will still be more readable in a BBCode code block than an attachment.

Or, if you have it online already in some accessible location (like a public github repository), you can post a link to that.

Also, if the code is really that long and complicated, expecting people to help you rewrite it in a forum discussion here might be too much to expect.

Ok understood. But I'm actually expecting them to copy and run it on their compiler and to provide me with the results.

 

[PeterDonis]

I'm actually expecting them to copy and run it on their compiler and to provide me with the results.

That might be too much to expect here too. You would need someone who was familiar not just with FORTRAN but with this particular code. The best place to find such a person would be at whatever place the code originally came from.

 

[member 657093]

Right, but something has to call these subroutines in the correct way to get answers. An example of that would really help. Someone might know FORTRAN very well (Not me. I am rusty and don't have a compiler.) and still not know how to use these subroutines.
If you know how to do it, even with your Python version, that would probably help a lot.

import numpy as nptt = 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)

# 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)# 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# time-dependent EMEP coefficients
M = 2.55e+19
O2 = 5.2e+18
XN2 = 1.99e+19
RPATH3 = 0.65
RPATH4 = 0.35# 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# 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)
# 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)

# 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)

 

[phyzguy]

Have you tried using f2py? It allows you to compile the fortran code as a Python module.
https://numpy.org/devdocs/f2py/f2py-examples.html

 

[Vanadium 50]

But I'm actually expecting them to copy and run it on their compiler and to provide me with the results.

Are you sure "expect" is the word you want? Demanding work of volunteers is seldom the best strategy, Other adjectives might be applied.

Besides:

I don't see anywhere that the final solution is printed.

So what you're really asking is for us to debug the FORTRAN code first, as it does not do what you want it to.

 

[member 657093]

Are you sure "expect" is the word you want? Demanding work of volunteers is seldom the best strategy, Other adjectives might be applied.

Besides:So what you're really asking is for us to debug the FORTRAN code first, as it does not do what you want it to.

I am actually asking help from a Fortran programmer to provide me with the values of the variables RC, DJ, Y, T, DY for TIME=331250 from the Fortran code as it stands to compare with my Python values to see if I understood the formulation correctly.

 

[Vanadium 50]

I think your best bet then is to hire a FORTRAN programmer. If the FORTRAN doesn't do exactly what you want it to do, it is not reasonable to expect us to modify it, much less demand it.

 

[jtbell]

I see several subroutines, but no main program. One of the subroutines (EMEPCF) is invoked (CALLed) by two of the others (FEVAL and JEVAL). None of the other subroutines (including FEVAL and JEVAL) are invoked anywhere in this code.

Clearly you need a main program that invokes these subroutines, in some specific sequence, and probably with some specific input data (unless some of these subroutines themselves initialize all the variables with data coded into them).

 

[Vanadium 50]

If the problem were to turn working FORTRAN into working Python, that would be one thing. But the problem is to turn non-working FORTRAN into working Python,

As people have noted, there is no main program, just subroutines. As the OP has noted, he needs more output than is provided by the code he posted.

Yes, I know he expects this. Doesn't mean its going to happen or is even possible.

 

[member 657093]

I think your best bet then is to hire a FORTRAN programmer. If the FORTRAN doesn't do exactly what you want it to do, it is not reasonable to expect us to modify it, much less demand it.

You are taking a reply I made to PeterDonis out of context! I didn't demand anything! I asked a question and I EXPECT someone who knows the subject to help. I didn't expect for others to play with words and make a whole issue out of the words!!!!!

 

[PeterDonis]

I EXPECT someone who knows the subject to help

Actually, "expect" in English does carry a connotation of "demand" as you are using it in the above quote.

As far as helping, sometimes the best help you will get is for people to point out to you that the problem is insoluble as you have described it. That might well be the case here, since the FORTRAN code isn't even complete to start with, so someone who is familiar with this particular FORTRAN library would have to complete the FORTRAN code first. The fact that no such person has come forward in this thread means there likely isn't such a person among PF's members. You would have to go back to the place where the FORTRAN code originally came from, as I suggested before.

 

[member 657093]

My problem is the following:

If I had a Python code that says:

def test(a):
    b = 13
    return a*b

The only addition I would need is to add:

c = test(5)

print(c)

As I am not Fortran literate, I was expecting some help in this way, where the Fortran programmer can call the subroutine and ask to print the values of the variables; just as I did for the second bit of Python code above.

 

[member 657093]

Actually, "expect" in English does carry a connotation of "demand" as you are using it in the above quote.

As far as helping, sometimes the best help you will get is for people to point out to you that the problem is insoluble as you have described it. That might well be the case here, since the FORTRAN code isn't even complete to start with, so someone who is familiar with this particular FORTRAN library would have to complete the FORTRAN code first. The fact that no such person has come forward in this thread means there likely isn't such a person among PF's members. You would have to go back to the place where the FORTRAN code originally came from, as I suggested before.

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!!

 

[PeterDonis]

My problem is the following:

If...

"If". But what is left out of the FORTRAN code actually that simple? What if it's not?

 

[member 657093]

"If". But what is left out of the FORTRAN code actually that simple? What if it's not?

I don't know that as I do not know Fortran! But are you saying that, as the code stands, a Fortran programmer cannot do anything to call up and compile a subroutine in that code?

 

[PeterDonis]

you took the second meaning

Yes, because in the context in which you used the term, the second meaning is the one a typical reader is more likely to infer. The first meaning is more likely in an impersonal context, for example: "I expect the sun to rise tomorrow morning." But if you say to me "I expect you to do X", I'm going to take that in the sense of the second meaning, i.e., that you think you have some reason or justification for asking me to do X for you, so that if I refuse, I have failed you in some way.

 

[PeterDonis]

I don't know that as I do not know Fortran!

Yes. And that means you don't know whether what you are asking for is actually that simple or not. What if it's not?

are you saying that, as the code stands, a Fortran programmer cannot do anything to call up and compile a subroutine in that code?

I am not a FORTRAN expert so I can't say. But people who apparently do know some FORTRAN are giving you responses in this thread that, to me, indicate that what you are asking for is not that simple and not something that someone unfamiliar with the FORTRAN library in question is going to be able to do.

 

[member 657093]

something that someone unfamiliar with the FORTRAN library in question is going to be able to do.

You told me to go to the source! But the source is only this! They have an introduction to the problem and what it involves. As the formulation is lengthy, they just pointed out to this Fortran code. source

 

[member 657093]

But people who apparently do know some FORTRAN

Who exactly? jtbell? Because to me he looks like he knows Fortran.

 

[member 657093]

you need a main program that invokes these subroutines

Can you please make one and call up these subroutines to print the values of variables RC, DJ, Y, T, DY for TIME=331250?

 

[FactChecker]

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.

 

[Mark44]

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.

      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.

 

[PeterDonis]

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.

 

[PeterDonis]

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.

 

[member 657093]

It would help if your Python indicated what it matches with a comment indicating the exact FORTRAN subroutine name.

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)
#-------------------------