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!!Filip Larsen said:Perhaps you can use https://www.onlinegdb.com/online_fortran_compiler or similar.
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 said: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.
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.Vick said: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.
You want my Python formulation?FactChecker said: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.
Please do not post code in attachments. You can use BBCode tags to put code directly in your post.Vick said:You want my Python formulation?
What if the code is very long?PeterDonis said:Please do not post code in attachments. You can use BBCode tags to put code directly in your post.
It will still be more readable in a BBCode code block than an attachment.Vick said:What if the code is very long?
Ok understood. But I'm actually expecting them to copy and run it on their compiler and to provide me with the results.PeterDonis said: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.
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.Vick said:I'm actually expecting them to copy and run it on their compiler and to provide me with the results.
FactChecker said: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)Are you sure "expect" is the word you want? Demanding work of volunteers is seldom the best strategy, Other adjectives might be applied.Vick said:But I'm actually expecting them to copy and run it on their compiler and to provide me with the results.
FactChecker said:I don't see anywhere that the final solution is printed.
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 said: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.
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!!!!!Vanadium 50 said: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.
Actually, "expect" in English does carry a connotation of "demand" as you are using it in the above quote.Vick said:I EXPECT someone who knows the subject to help
def test(a):
b = 13
return a*bc = test(5)
print(c)PeterDonis said: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.
"If". But what is left out of the FORTRAN code actually that simple? What if it's not?Vick said:My problem is the following:
If...
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 said:"If". But what is left out of the FORTRAN code actually that simple? What if it's not?
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.Vick said:you took the second meaning
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?Vick said:I don't know that as I do not know Fortran!
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.Vick said:are you saying that, as the code stands, a Fortran programmer cannot do anything to call up and compile a subroutine in that code?
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. sourcePeterDonis said:something that someone unfamiliar with the FORTRAN library in question is going to be able to do.
Who exactly? jtbell? Because to me he looks like he knows Fortran.PeterDonis said:But people who apparently do know some FORTRAN
Can you please make one and call up these subroutines to print the values of variables RC, DJ, Y, T, DY for TIME=331250?jtbell said:you need a main program that invokes these subroutines
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.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!!
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.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?
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)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.Vick said:You told me to go to the 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.Vick said:As the formulation is lengthy, they just pointed out to this Fortran code. source
FactChecker said: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)
#-------------------------