-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathCondL_h.m
51 lines (43 loc) · 2.06 KB
/
CondL_h.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
function CondL_h
global Theta_LL Theta_r Theta_s Alpha hh n m Se KL_h Ks MN ML ND DTheta_LLh
global NL J Theta_L h IS KIT TT Thmrlefc CKT CKTN POR
% PRN The lowest suction head (The maximum value of matric head,considering
% the negative sign before them. The absolute value of which is smallest) at which soil remains saturated.
% PRN=-1e-6;
MN=0;
for ML=1:NL
J=IS(ML);
for ND=1:2
MN=ML+ND-1;
if hh(MN)>=-1e-6
Theta_LL(ML,ND)=Theta_s(J);
hh(MN)=-1e-6;
DTheta_LLh(ML,ND)=0;
Se(ML,ND)=1;
elseif hh(MN)<=-1e5
Theta_LL(ML,ND)=Theta_r(J);
hh(MN)=-1e5;
DTheta_LLh(ML,ND)=0;
Se(ML,ND)=0;
else
Theta_LL(ML,ND)=Theta_r(J)+(Theta_s(J)-Theta_r(J))/(1+abs(Alpha(J)*hh(MN))^n(J))^m(J);
if Thmrlefc
DTheta_LLh(ML,ND)=(Theta_s(J)-Theta_r(J))*Alpha(J)*n(J)*abs(Alpha(J)*hh(MN))^(n(J)-1)*(-m(J))*(1+abs(Alpha(J)*hh(MN))^n(J))^(-m(J)-1);
else
if abs(hh(MN)-h(MN))<1e-3
DTheta_LLh(ML,ND)=(Theta_s(J)-Theta_r(J))*Alpha(J)*n(J)*abs(Alpha(J)*hh(MN))^(n(J)-1)*(-m(J))*(1+abs(Alpha(J)*hh(MN))^n(J))^(-m(J)-1);
else
DTheta_LLh(ML,ND)=(Theta_LL(ML,ND)-Theta_L(ML,ND))/(hh(MN)-h(MN));
end
end
Se(ML,ND)=Theta_LL(ML,ND)/POR(J);
end
if KIT
CKT(MN)=CKTN/(50+2.575*TT(MN));
KL_h(ML,ND)=CKT(MN)*Ks(J)*(Se(ML,ND)^(0.5))*(1-(1-Se(ML,ND)^(1/m(J)))^m(J))^2;
end
end
end
%%%%%%%%% Unit of KL_h is determined by Ks, which would be given at the%%%%
%%%%%%%%% beginning.Import thing is to keep the unit of matric head hh(MN)
%%%%%%%%% as 'cm'.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%