-
Notifications
You must be signed in to change notification settings - Fork 0
/
modelinputs_v12_1b.bak
232 lines (200 loc) · 8.62 KB
/
modelinputs_v12_1b.bak
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
229
230
231
% === MyLake model, version 1.2, 15.03.05 === %%%% KOJI very quick 221110
% by Tom Andersen & Tuomo Saloranta, NIVA 2005
%
% Module for reading input data and parameters
% Code checked by TSA, xx.03.2005
% Last modified by TSA, 15.08.2006 (Az replaced by In_Az 10.03.06; Possibility to have NaN in Global rad. series, 15.08.06)
function [In_Z,In_Az,tt,In_Tz,In_Cz,In_Sz,In_TPz,In_DOPz,In_Chlz,In_DOCz,In_TPz_sed,In_Chlz_sed,In_O2z,In_NO3z,In_NH4z,In_SO4z,In_Fe2z,In_Ca2z,In_pHz,In_CH4z,In_Fe3z,In_Al3z,In_SiO4z,In_SiO2z,In_diatomz,In_FIM,Ice0,Wt,Inflw,...
Phys_par,Phys_par_range,Phys_par_names,Bio_par,Bio_par_range,Bio_par_names] ...
= modelinputs_v12_1b(M_start,M_stop,init_filename,init_sheet,...
input_filename,input_sheet,param_filename,param_sheet,dt);
% Inputs:
% M_start : Model start date [year, month, day]
% M_stop : Model stop date [year, month, day]
% + Input filenames and sheetnames
% dt : Time step (= 1 day)
% Outputs:
% tt : Solution time domain (day)
% In_Z : Depths read from initial profiles file (m)
% In_Az : Areas read from initial profiles file (m2)
% In_Tz : Initial temperature profile read from initial profiles file (deg C)
% In_Cz : Initial tracer profile read from initial profiles file (-)
% In_Sz : Initial sedimenting tracer (or suspended inorganic matter) profile read from initial profiles file (kg m-3)
% In_TPz : Initial total P profile read from initial profiles file (mg m-3)
% In_DOPz : Initial dissolved organic P profile read from initial profiles file (mg m-3)
% In_Chlz : Initial chlorophyll a profile read from initial profiles file (mg m-3)
% In_DOCz : Initial DOC profile read from initial profiles file (mg m-3)
% In_TPz_sed : Initial total P profile in the sediment compartments read from initial profiles file (mg m-3)
% In_Chlz_sed : Initial chlorophyll a profile in the sediment compartments read from initial profiles file (mg m-3)
% ...
% In_FIM : Initial profile of volume fraction of inorganic matter in the sediment solids (dry weight basis)
% Ice0 : Initial conditions, ice and snow thicknesses (m) (Ice, Snow)
% Wt : Weather data
% Inflow : Inflow data
% Phys_par : Main 23 parameters that are more or less fixed
% Phys_par_range : Minimum and maximum values for Phys_par (23 * 2)
% Phys_par_names : Names for Phys_par
% Bio_par : Main 15 parameters that are more or less site specific
% Bio_par_range : Minimum and maximum values for Bio_par (15 * 2)
% Bio_par_names : Names for Bio_par
global ies80;
% == Read model parameter file
% ParaMx = xlsread(param_filename);
% % Koji 06.10.10
% ParaMx = dlmread(param_filename, '\t', 'B3..D48'); % 46 by 3 This is
% % version specific
f = fopen(param_filename);
garbage = fgetl(f);
garbage = fgetl(f); % file get line
data = textscan(f, '%s%f%f%f%s', 60, 'Delimiter', '\t');
fclose(f);
% read 46 lines
par_names = data{1};
par = data{2};
par_range = [data{3} data{4}]; % horizontally concatenate
Phys_par_names = par_names(1:23,:);
Bio_par_names = par_names(24:60,:);
Phys_par = par(1:23,:);
Bio_par = par(24:60,:);
Phys_par_range = par_range(1:23,:);
Bio_par_range = par_range(24:60,:);
% Bio_par(34) = 0.25;
% Bio_par(35) = 1.95;
% Bio_par(36) = 2.65;
% %Main physical parameters (dz, Kz_ak, etc...)
% % Phys_par_names=ParaMx(3:25,1);
% % Phys_par=ParaMx(3:25,2);
% % Phys_par_range=ParaMx(3:25,3:4);
% Phys_par_names=NaN;
% Phys_par=ParaMx(1:23,1);
% Phys_par_range=ParaMx(1:23,2:3);
% %Main biological parameters (Y_cp, m_twty, g_twty, etc...)
% % Bio_par_names=ParaMx(26:48,1);
% % Bio_par=ParaMx(26:48,2);
% % Bio_par_range=ParaMx(26:48,3:4);
% Bio_par_names=NaN;
% Bio_par=ParaMx(24:46,1);
% Bio_par_range=ParaMx(24:46,2:3);
% Vertical settling velocities
U = Bio_par(8:9); %for sedimenting velocities
if any(U<0)
error('Given settling velocity must be positive')
end
% == Read morphometric and initial profile file
% InitMx = xlsread(init_filename);
% % Koji 06.10.10
InitMx = dlmread(init_filename, '\t', 2, 0);
%% changed below 3:end to 1:end
In_Z=InitMx(1:end,1);
In_Az=InitMx(1:end,2);
In_Tz=InitMx(1:end,3);
In_Cz=InitMx(1:end,4);
In_Sz=InitMx(1:end,5);
In_TPz=InitMx(1:end,6);
In_DOPz=InitMx(1:end,7);
In_Chlz=InitMx(1:end,8);
In_DOCz=InitMx(1:end,9);
In_TPz_sed=InitMx(1:end,10);
In_Chlz_sed=InitMx(1:end,11);
In_FIM=InitMx(1:end,12);
Ice0=InitMx(1,13:14);
In_O2z=InitMx(1:end,15);
In_NO3z=InitMx(1:end,16);
In_NH4z=InitMx(1:end,17);
In_SO4z=InitMx(1:end,18);
In_HSz=InitMx(1:end,19);
In_H2Sz=InitMx(1:end,20);
In_Fe2z=InitMx(1:end,21);
In_Ca2z=InitMx(1:end,22);
In_pHz=InitMx(1:end,23);
In_CH4z=InitMx(1:end,24);
In_Fe3z=InitMx(1:end,25);
In_Al3z=InitMx(1:end,26);
In_SiO4z=InitMx(1:end,27);
In_SiO2z=InitMx(1:end,28);
In_diatomz=InitMx(1:end,29);
tt = [datenum(M_start):dt:datenum(M_stop)]'; % Solution time domain
% == Read input forcing data file
% $$$ % KOJI
% $$$ if findstr(input_filename, 'csv')
% $$$ InputMx = csvread(input_filename);
% $$$ SIZEofINPUT = size(InputMx);
% $$$ cushion = zeros(2, SIZEofINPUT(1, 2));
% $$$ InputMx = [cushion; InputMx];
% $$$ else
% $$$ [InputMx,StrMx]=xlsread(input_filename,input_sheet);
% $$$ end
% InputMx = xlsread(input_filename);
% % Koji 06.10.10
InputMx = dlmread(input_filename, '\t', 2, 0);
%% changed below 3:end to 1:end 2011-09-29
In_Date=InputMx(1:end,1:3);
In_Met=InputMx(1:end,4:10);
In_Inflow=InputMx(1:end,11:32);
tmet=datenum(In_Date);
dum=100*((tmet(end)-tmet(1)+1)-length(tmet))/(tmet(end)-tmet(1)+1);
% $$$ disp(['Percent missing dates in meteorology and inflow data: ']);
% $$$ disp([num2str(dum) ' %']);
dum=100*sum(isnan(In_Met))./length(tmet);
% $$$ disp(['Percent missing values in meteorology data (values correspond to columns 4-10 in input file): ']);
% $$$ disp([num2str(dum) ' %']);
dum=100*sum(isnan(In_Inflow))./length(tmet);
% $$$ disp(['Percent missing values in inflow data (values correspond to columns 11-17 in input file): ']);
% $$$ disp([num2str(dum) ' %']);
% $$$ disp(' ')
clear Wt
for i=1:7 %Interpolate over missing values and dates
nonnans = find(isnan(In_Met(:,i))==0);
if(isempty(nonnans)) % if the whole column is NaNs then preserve it
Wt(:,i) = NaN*ones(length(tt(:)),1);
else
repaired = interp1(nonnans,In_Met(nonnans,i),[1:length(In_Met(:,i))]);
Wt(:,i) = interp1(tmet, repaired, tt(:));
end
end
% Wt(:,1) Global radiation (MJ/(m^2 day))
% Wt(:,2) Cloud cover (-)
% Wt(:,3) Air temperature (deg. C, at 2 m height)
% Wt(:,4) Relative humidity (%, at 2 m height)
% Wt(:,5) Air pressure (mbar)
% Wt(:,6) Wind speed (m/s at 10 m height)
% Wt(:,7) Precipitation (mm/day)
clear Inflw
for i=1:22 %Interpolate over missing values and dates
nonnans = find(isnan(In_Inflow(:,i))==0);
if(isempty(nonnans)) % if the whole column is NaNs then preserve it
Inflw(:,i) = NaN*ones(length(tt(:)),1);
else
repaired = interp1(nonnans,In_Inflow(nonnans,i),[1:length(In_Inflow(:,i))]);
Inflw(:,i) = interp1(tmet, repaired, tt(:));
end
end
%Inflw(:,9) = Inflw(:,8); %DIC
%Inflw(:,10) = (Inflw(:,8)+1) ./ (Inflw(:,8)+1) .* 8000 ; % O2 Must be defined better
% Inflw(:,1) Inflow volume (m3 day-1)
% Inflw(:,2) Inflow temperature (deg C)
% Inflw(:,3) Inflow tracer concentration (-)
% Inflw(:,4) Inflow sedimenting tracer (or suspended inorganic matter) concentration (kg m-3)
% Inflw(:,5) Inflow total phosphorus (TP) concentration (mg m-3)
% Inflw(:,6) Inflow dissolved organic phosphorus (DOP) concentration (mg m-3)
% Inflw(:,7) Inflow chlorophyll a concentration (mg m-3)
% Inflw(:,8) Inflow DOC concentration (mg m-3)
% This will be created after fact ... RMC 10.03.2014
% Inflow DIC concentration (mg m-3)
% Inflw(:,10) Inflow O2 concentration (mg m-3)
% International Equation of State 1980
% 5-order polynomial for density as function of temperature
ies80 = [6.536332e-9,-1.120083e-6,1.001685e-4,-9.09529e-3,6.793952e-2,999.842594];
% Default turbulence and wind shelter parameterization (Hondzo and Stefan, 1993; Ellis et al., 1991)
if(isnan(Phys_par(2)))
Phys_par(2) = 0.00706*(In_Az(1)/1e6)^0.56; % default diffusion coeff. parameterisation
end
if(isnan(Phys_par(3)))
Phys_par(3) = 8.98e-4; %default value for diffusion coeff. in ice-covered water
end
if(isnan(Phys_par(4)))
Phys_par(4) = 7e-5; % default minimum allowed stability frequency, N2 > N0 <=> Kz < Kmax (1/s2)
end
if(isnan(Phys_par(5)))
Phys_par(5) = 1-exp(-0.3*In_Az(1)/1e6); % default wind sheltering parameterisation
end