All variables, constants, parameters, etc.

double KPF,KNF,KFZ,KDZ,KFB,KZB,KDB,KZC,KBC,KAC,KFS,KDS, MUFMAX,MUFZMAX,MUDZMAX,MUFBMAX,MUZBMAX, MUDBMAX,MUZCMAX,MUBCMAX,MUACMAX,MUFSMAX,MUDSMAX, MBF,MBZ,MBB,MBC,MBS,MBBC,MBBS, RCMAX,RSMAX,MBOF,MBOZ,MBOB,MBOC,MBOS, MF,MZ,MB,RESPF, RESPZ,RESPB,RESPC,RESPS, UDP,UDN,PHOT,ALPHA,RE,OKT,SED,SEDP,LOPT, KW,KF,KD,KPD,HK,BKR,ZCR,FCR,MB2,MZ2,MF2, LABDAB,LABDAZ,LABDAF, KA,MK,COP,MG,COD, HLP1,HLP2,HLP3,HLP4,HLP5,HLP6,HLP7,HLP8, HLP9,HLP10,HLP11,HLP12, HLP13,HLP14,HLP15,HLP16, QGF,QFZ,QFB,QFS,QFD, QFE,QDZ,QZB, QZC,QZD,QZE,QDB,QBC, QBD,QBE,QAC,QDS,QCD,QCE, QSD,QSE,QDP, QPF,QPD,QDN,QNF, QFO,QOF,QOZ, QOB,QOC, QOS,QOD,QAD, RC,RS,L,F,D,Z,S,B,C,A,P,N,F0,HKTEMP, DELTAPU,DELTANU, DELTAOU,DELTAAU, PU,NU,OU,AU, DFDT,DZDT,DBDT,DCDT, DSDT,DPDT,DNDT, DODT,DADT,DDDT, PTOEG,NTOEG,OTOEG,VTOEG, O,REA, TEMP,L_JAARMAX,TEMP_JAARMAX, Error;
PKWAD,NKWAD,DPU,DNU, t,dt,Tmin,Tmax.

All starting values

KPF = 0.1;
KNF = 0.1;
KFZ = 15.0;
KDZ = 60.0;
KFB = 15.0;
KZB = 1.0;
KDB = 60.0;
KZC = 1.0;
KBC = 5.0;
KAC = 0.2;
KFS = 20.0;
KDS = 60.0;

MUFMAX = 4.5;
MUFZMAX = 1.4;
MUDZMAX = 0.5;
MUFBMAX = 0.2;
MUZBMAX = 0.4;
MUDBMAX = 0.2;
MUZCMAX = 0.02;
MUBCMAX = 0.06;
MUACMAX = 0.03;
MUFSMAX = 0.1;
MUDSMAX = 0.07;

MBF = 0.3;
MBZ = 0.3;
MBB = 0.3;
MBC = 0.3;
MBS = 0.3;

MBBC = 0.4;
MBBS = 0.4;

RCMAX = 13.0;
RSMAX = 10.0;

MBOF = 0.001;
MBOZ = 0.001;
MBOB = 0.001;
MBOC = 0.001;
MBOS = 0.001;

MF = 0.025;
MZ = 0.005;
MB = 0.05;

RESPF = 0.001;
RESPZ = 0.01;
RESPB = 0.01;
RESPC = 0.01;
RESPS = 0.01;

UDP = 0.00004;
UDN = 0.002;
PHOT = 1.0;
ALPHA = 0.2;
RE = 0.4;
OKT = 0.085;

SED = 0.05;
SEDP = 0.1;
LOPT = 3000.0;

KW = 0.2;
KF = 0.05;
KD = 0.8;

KPD = 0.5;
HK = 0.1;
BKR = 20.0;
ZCR = 5.0;
FCR = 30.0;
MB2 = 10.0;
MZ2 = 3.0;
MF2 = 15.0;

LABDAB = 1.0;
LABDAZ = 1.0;
LABDAF = 1.0;

KA = 2.0;
MK = 5.0;
COP = 2.0;
MG = 1.0;
COD = 2.0;

HLP1 = MBBC / RCMAX;
HLP2 = MBBS / RSMAX;
HLP3 = XKWAD (KPF);
HLP4 = XKWAD (KNF);
HLP5 = XKWAD (KFZ);
HLP6 = XKWAD (KFB);
HLP7 = XKWAD (KFS);
HLP8 = XKWAD (KDZ);
HLP9 = XKWAD (KZB);
HLP10 = XKWAD (KBC);
HLP11 = W(MUBCMAX, HLP10, BKR);
HLP12 = XKWAD (KZC);
HLP13 = XKWAD (KDB);
HLP14 = XKWAD (KAC);
HLP15 = XKWAD (KDS);

Topt[1]=24.0;
Topt[2]=24.0;
Topt[3]=24.0;
Topt[4]=26.0;
Topt[5]=26.0;
Topt[6]=26.0;

Q1[1]=15.0;//=24-9
Q1[2]=15.0;
Q1[3]=15.0;
Q1[4]=14.0;
Q1[5]=12.5;
Q1[6]=11.0;

Q2[1]=10.0;//=34-24
Q2[2]=10.0;
Q2[3]=10.0;
Q2[4]=11.0;
Q2[5]=11.0;
Q2[6]=11.0;

M[1]=3.0;//dummy
M[2]=3.0;//zuurstofhalfwaarde voor onderhoud Z
M[3]=3.0;//dummy
M[4]=3.0;//zuurstof ... C
M[5]=3.0;//zuurstof ... S
M[6]=3.0;//zuurstof ... H

LABDA[1]=1.0;//dummy
LABDA[2]=1.0;//helling zuurstofcurve Z
LABDA[3]=1.0;//dummy
LABDA[4]=1.0;//helling zuurstofcurve C
LABDA[5]=1.0;//helling zuurstofcurve S
LABDA[6]=1.0;//helling zuurstofcurve H

F = 20.0;
Z = 0.0;
B = 0.0;
C = 0.0;
S = 0.0;
P = 0.03;
N = 0.15;
O = 9.8;
A = 0.0;
D = 7.2;
L = 0.0;

F0 = 2.0;

L_JAARMAX = 3000.0;
TEMP_JAARMAX = 20.0;

REA = REA_MOMENT();
TEMP = TEMP_MOMENT();

NU = 0.0;
PU = 0.0;
OU = 0.0;
AU = 0.0;

DPU = 0.01;
DNU = 0.03;

PTOEG = 0.0;
NTOEG = 0.0;
OTOEG = 0.0;
VTOEG = 0.0;
Tmin = 0.0;
T = 0.0;
dt = 0.04;
Tmax = 30.1;
}

Three input variables

double TEMP_MOMENT()
{
double TEMP_TIJD, TEMP_SEIZOEN, resultaat;
TEMP_TIJD = -Math.cos(6.283185 * (T - Math.floor(T)));
TEMP_SEIZOEN = 6+6+((1-Math.cos(t*0.0172142))*(TEMP_JAARMAX-7)*0.5);
resultaat = TEMP_SEIZOEN * TEMP_TIJD;
return resultaat;
}

double L_MOMENT()
{
double L_NU, L_TIJD, L_SEIZOEN, resultaat;

L_TIJD = - Math.cos( 6.283185 * (T - Math.floor(T)));
L_SEIZOEN = - Math.cos(T * 0.0172142) * 0.5;
L_NU = (L_TIJD + L_SEIZOEN) * L_JAARMAX / 1.5;
if (L_NU < 0.0) {resultaat = 0.0;}
else {resultaat = L_NU;}
return resultaat;
}

double REA_MOMENT()
{
double TEMPTEMP, TEMPTEMPTEMP, resultaat;
TEMPTEMP = XKWAD(TEMP);
TEMPTEMPTEMP = TEMPTEMP * TEMP;
resultaat = 14.61996 -0.40420 * TEMP + 0.00842 * TEMPTEMP - 0.00009 * TEMPTEMPTEMP;
return resultaat;
}

The model in the loop

Do ... i to ...

PKWAD = XKWAD(P);
NKWAD = XKWAD(N);

QGF = FT(1)*FF(L,F,D)*MUFMAX*MIN(PKWAD/(HLP3+PKWAD),NKWAD/(5*(HLP4 +NKWAD)))*F;
QFZ = FT(2)*FO(2) * W(MUFZMAX,HLP5,F)*Z;
QFB = FT(3)*W(MUFBMAX,HLP6,F)*B;
QFS = FT(5)*FO(5) * S*W(MUFSMAX,HLP7,F);
QFD = MBF*QGF;
if (F > F0) {QFD = QFD + MF*F;};
QFE = MBOF*F;

QDZ = FT(2)*FO(2) * Z*W(MUDZMAX,HLP8,D);
QZB = FT(3) * W(MUZBMAX,HLP9,Z)*B;
HLP16 = W(MUBCMAX,HLP10,B);

QZC = FT(4)*FO(4)*MIN((HLP11-HLP16),(W(MUZCMAX,HLP12,Z)*JOTA(B,LABDAB,MB2)))*C;
QZD = MBZ*(QFZ + QDZ) + FOX(O)*MZ*Z;
QZE = MBOZ*Z;
QDB = FT(3)*B*W(MUDBMAX,HLP13,D);
QBC = FT(4)*FO(4)*HLP16*C;
QBD = MBB * (QFB + QZB + QDB);
QBE = MBOB*B;
QAC = FT(4) * FO(4) * C*W(MUACMAX,HLP14,A);
QDS = FT(5) * FO(5) * S*W(MUDSMAX,HLP15,D);
RC = QAC * QBC + QZC;
RS = QFS + QDS;
QCD = (MBC + HLP1 * RC) * RC;
QCE = MBOC * C;
QSD = (MBS + HLP2 * RS) * RS;
QSE = MBOS * S;
QDP = UDP * E1(TEMP) * E2(O) * D;
QPF = QGF / 106.0;
QPD = SEDP * E3(O) * P;
QDN = UDN * E1(TEMP) * E2(O) * D;
QNF = QGF * 0.0472;
QFO = PHOT * QGF;
QOF = RESPF * F;
QOZ = RESPZ * Z;
QOB = RESPB * B;
QOC = RESPC * C;
QOS = RESPS * S;
QOD = OKT * D;
QAD = ALPHA * A;

DELTAPU = PU * dt;
DELTANU = NU * dt;
DELTAOU = OU * dt;
DELTAAU = AU * dt;

DFDT = (QGF - QFZ - QFB - QFS - QFD - QFE)*dt;
DZDT = (QFZ + QDZ - QZB - QZC - QZD - QZE)*dt;
DBDT = (QFB + QZB + QDB - QBC - QBD - QBE)*dt;
DCDT = (QAC + QBC + QZC - QCD - QCE)*dt;
DSDT = (QFS + QDS - QSD - QSE)*dt;

DPDT = DELTAPU + (QDP - QPF - QPD + DPU)*dt;
DNDT = DELTANU + (QDN - QNF + DNU)*dt;
DODT = DELTAOU + (QFO + RE * (REA - O) - QOF - QOZ - QOB - QOC - QOS - QOD)*dt;
DADT = DELTAAU + (AU - QAC - QAD)*dt;
DDDT = (QFD +QZD +QBD +QCD +QSD +QAD -QPD -QDN -QDZ -QDB -QDS -SED*D)*dt;

F = F + DFDT; if (F < 0.0) {F = 0.0;}
Z = Z + DZDT; if (Z < 0.0) {Z = 0.0;}
B = B + DBDT; if (B < 0.0) {B = 0.0;}
C = C + DCDT; if (C < 0.0) {C = 0.0;}
S = S + DSDT; if (S < 0.0) {S = 0.0;}
P = P + DPDT; if (P < 0.0) {P = 0.0;}
N = N + DNDT; if (N < 0.0) {N = 0.0;}
O = O + DODT; if (O < 0.0) {O = 0.0;}
A = A + DADT; if (A < 0.0) {A = 0.0;}
D = D + DDDT; if (D < 0.0) {D = 0.0;}
PTOEG = PTOEG + DELTAPU;
NTOEG = NTOEG * DELTANU;
OTOEG = OTOEG + DELTAOU;
VTOEG = VTOEG + DELTAAU;

Graph1.newValue(N);
Graph2.newValue(TEMP);
Graph3.newValue(QDP);
Graph4.newValue(E1(TEMP));
Graph5.newValue(O);
Graph6.newValue(F);
Graph7.newValue(P);
Graph8.newValue(E2(O));

Text1.setText(Double.toString(N));
Text2.setText(Double.toString(TEMP));
Text3.setText(Double.toString(QDP));
Text4.setText(Double.toString(E1(TEMP)));
Text5.setText(Double.toString(O));
Text6.setText(Double.toString(F));
Text7.setText(Double.toString(P));
Text8.setText(Double.toString(E2(O)));

until (T < Tmax);

Enschede, 16/dec/2001