forked from yijianzeng/STEMMUS
-
Notifications
You must be signed in to change notification settings - Fork 0
/
MainLoop.m
228 lines (215 loc) · 6.61 KB
/
MainLoop.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
clc;
clearvars -except TRAP_ind EVAP_ind Sim_Temp_ind Sim_Theta_ind;
close all;
run Constants
if Evaptranp_Cal==1
run EvapTransp_Cal;
end
% function MainLoop
global i tS KT Delt_t TEND TIME MN NN NL ML ND hOLD TOLD h hh T TT P_gOLD P_g P_gg Delt_t0
global KIT NIT TimeStep Processing
global SUMTIME hhh TTT P_ggg Theta_LLL DSTOR Thmrlefc CHK Theta_LL Theta_L
global NBCh AVAIL Evap DSTOR0 EXCESS QMT RS BCh hN hSAVE NBChh DSTMAX Soilairefc Trap sumTRAP_dir sumEVAP_dir
global TSAVE IRPT1 IRPT2 AVAIL0 TIMEOLD TIMELAST SRT ALPHA BX alpha_h bx Srt
global QL QL_h QL_T QV Qa KL_h Chh ChT Khh KhT TRAP_ind
run StartInit; % Initialize Temperature, Matric potential and soil air pressure.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TIMEOLD=0;
TIMELAST=0;
for i=1:tS+1 % Notice here: In this code, the 'i' is strictly used for timestep loop and the arraies index of meteorological forcing data.
KT=KT+1 % Counting Number of timesteps
if KT>1 && Delt_t>(TEND-TIME)
Delt_t=TEND-TIME; % If Delt_t is changed due to excessive change of state variables, the judgement of the last time step is excuted.
end
TIME=TIME+Delt_t; % The time elapsed since start of simulation
TimeStep(KT,1)=Delt_t;
SUMTIME(KT,1)=TIME;
Processing=TIME/TEND
%%%%%% Updating the state variables. %%%%%%%%%%%%%%%%%%%%%%%%%%%%
if TIME>=489*3600 && TIME<=492*3600 %7-13 9-10 p=52mm
NBChh=1;
elseif TIME>=952*3600 && TIME<=964*3600 %8-1 16-18 p=60mm
NBChh=1;
elseif TIME>=1288*3600 && TIME<=1297*3600 %8-15 16-17 p=67mm
NBChh=1;
elseif TIME>=1862*3600 && TIME<=1870*3600 %9-8 14-18 p=93.11mm
NBChh=1;
else
NBChh=2;
end
if IRPT1==0 && IRPT2==0
for MN=1:NN
hOLD(MN)=h(MN);
h(MN)=hh(MN);
hhh(MN,KT)=hh(MN);
KL_h(MN,KT)=KL_h(MN,2);
Chh(MN,KT)=Chh(MN,2);
ChT(MN,KT)=ChT(MN,2);
Khh(MN,KT)=Khh(MN,2);
KhT(MN,KT)=KhT(MN,2);
if Thmrlefc==1
TOLD(MN)=T(MN);
T(MN)=TT(MN);
TTT(MN,KT)=TT(MN);
end
if Soilairefc==1
P_gOLD(MN)=P_g(MN);
P_g(MN)=P_gg(MN);
P_ggg(MN,KT)=P_gg(MN);
end
if rwuef==1
SRT(MN,KT)=Srt(MN,1);
ALPHA(MN,KT)=alpha_h(MN,1);
BX(MN,KT)=bx(MN,1);
end
end
DSTOR0=DSTOR;
if KT>1
run SOIL1
end
end
if Delt_t~=Delt_t0
for MN=1:NN
hh(MN)=h(MN)+(h(MN)-hOLD(MN))*Delt_t/Delt_t0;
TT(MN)=T(MN)+(T(MN)-TOLD(MN))*Delt_t/Delt_t0;
end
end
hSAVE=hh(NN);
TSAVE=TT(NN);
if NBCh==1
hN=BCh;
hh(NN)=hN;
hSAVE=hN;
elseif NBCh==2
if NBChh~=2
if BCh<0
hN=DSTOR0;
hh(NN)=hN;
hSAVE=hN;
else
hN=-1e6;
hh(NN)=hN;
hSAVE=hN;
end
end
else
if NBChh~=2
hN=DSTOR0;
hh(NN)=hN;
hSAVE=hN;
end
end
run Forcing_PARM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for KIT=1:NIT % Start the iteration procedure in a time step.
run SOIL2;
run CondL_T;
run Density_V;
run CondL_Tdisp;
run Latent;
run Density_DA;
run CondT_coeff;
run Condg_k_g;
run CondV_DE;
run CondV_DVg;
run h_sub;
if NBCh==1
DSTOR=0;
RS=0;
elseif NBCh==2
AVAIL=-BCh;
EXCESS=(AVAIL+QMT(KT))*Delt_t;
if abs(EXCESS/Delt_t)<=1e-10,EXCESS=0;end
DSTOR=min(EXCESS,DSTMAX);
RS=(EXCESS-DSTOR)/Delt_t;
else
AVAIL=AVAIL0-Evap(KT);
EXCESS=(AVAIL+QMT(KT))*Delt_t;
if abs(EXCESS/Delt_t)<=1e-10,EXCESS=0;end
DSTOR=0;
RS=0;
end
if Soilairefc==1
run Air_sub;
end
if Thmrlefc==1
run Enrgy_sub;
end
if max(CHK)<0.001
break
end
hSAVE=hh(NN);
TSAVE=TT(NN);
end
TIMEOLD=KT;
KIT
KIT=0;
run SOIL2;
% run TimestepCHK
if IRPT1==0 && IRPT2==0
if KT % In case last time step is not convergent and needs to be repeated.
MN=0;
for ML=1:NL
for ND=1:2
MN=ML+ND-1;
Theta_LLL(ML,ND,KT)=Theta_LL(ML,ND);
Theta_L(ML,ND)=Theta_LL(ML,ND);
end
end
run ObservationPoints
end
if (TEND-TIME)<1E-3
for MN=1:NN
hOLD(MN)=h(MN);
h(MN)=hh(MN);
hhh(MN,KT)=hh(MN);
if Thmrlefc==1
TOLD(MN)=T(MN);
T(MN)=TT(MN);
TTT(MN,KT)=TT(MN);
end
if Soilairefc==1
P_gOLD(MN)=P_g(MN);
P_g(MN)=P_gg(MN);
P_ggg(MN,KT)=P_gg(MN);
end
end
break
end
end
for MN=1:NN
QL(MN,KT)=QL(MN);
QL_h(MN,KT)=QL_h(MN);
QL_T(MN,KT)=QL_T(MN);
Qa(MN,KT)=Qa(MN);
QV(MN,KT)=QV(MN);
end
end
% run PlotResults
if Evaptranp_Cal==1 % save the variables for ETind scenario
Sim_Theta_ind=Sim_Theta;
Sim_Temp_ind=Sim_Temp;
TRAP=36000.*Trap;
TRAP_ind=TRAP';
EVAP=36000.*Evap;
EVAP_ind=EVAP';
disp ('Convergence Achieved for ETind scenario. Please switch to ETdir scenario and run again.')
else
TRAP=36000.*Trap;
TRAP_dir=TRAP';
EVAP=36000.*Evap;
EVAP_dir=EVAP';
for i=1:KT/24
% sumTRAP_ind(i)=0;
% sumEVAP_ind(i)=0;
sumTRAP_dir(i)=0;
sumEVAP_dir(i)=0;
for j=(i-1)*24+1:i*24
% sumTRAP_ind(i)=TRAP_ind(j)+sumTRAP_ind(i);
% sumEVAP_ind(i)=EVAP_ind(j)+sumEVAP_ind(i);
sumTRAP_dir(i)=TRAP(j)+sumTRAP_dir(i);
sumEVAP_dir(i)=EVAP(j)+sumEVAP_dir(i);
end
end
run PlotResults1
end