- Set up a run folder for the job (i.e., test_run01), and put executable file (i.e., fvcom), namelist file (superior_run.nml, only_DA.nml), set up an input folder for input files (i.e., superior_xxx.dat files and the atmospheric forcing file), and also set up an output folder.
- In only_DA.nml, turn on two sections:
(1) For lake surface temperature
(2) For vertical temperature
- In input folder, it should contain:
- all_obs_ts.xy
3 # the number of total mooring points
01 817036 4637430.5 8.00 08 0.0000 0.0000
# mooring number; mooring longitude; mooring latitude; water depth; layer number; 0; 0
1.7000
2.7000
4.7000
5.7000
6.2000
6.7000
7.2000
7.7000
02 871864.75 4637383 11.50 02 0.0000 0.0000
3.5000
10.5000
05 944111.5 4654256.5 24.50 17 0.0000 0.0000
1.0000
3.0000
5.0000
7.0000
9.0000
11.0000
13.0000
15.0000
16.5000
17.0000
17.5000
18.5000
19.5000
20.5000
21.5000
22.5000
23.5000
The matlab code used to generate “all_obs_ts.xy”
clc;clear;close;
% load the layer depth of observation data
ncload('/Users/chuyanzhao/Desktop/gls_observation_data/obs_glerl_mooring/glerl_lake_superior_temperature_mooring_2018_2020.nc','z','lat','lon');
% define variables
obs_num=1;
depth=196.6; % water depth at mooring location
layer_num=17; % layers number used in TS DA
filename=('all_obs_ts.xy');
fileID = fopen(filename,'w');
for i=1:obs_num
obs_pos(i,1)=lon+360;
obs_pos(i,2)=lat;
end
% print
for i=1:obs_num
fprintf(fileID,'%d\n',i);
fprintf(fileID,'%d %8.4f %8.4f %6.3f %d %8.4f %8.4f\n',i,obs_pos(i,1),obs_pos(i,2),depth,layer_num,0,0);
fprintf(fileID,'%8.4f\n',7,z(6:21));
end
fclose(fileID);
- all_obs_ts.dat
# mooring number; times
1 1866
# time; temperature; salnity
53544.8008 21.0049992 0.00000000 20.7539997 0.00000000 21.0209999 0.00000000 53544.8242 21.0049992 0.00000000 20.7539997 0.00000000 21.0209999 0.00000000 53544.8438 21.0680008 0.00000000 20.8799992 0.00000000 21.0839996 0.00000000
53544.8633 21.0680008 0.00000000 20.7539997 0.00000000 21.0839996 0.00000000 53544.8867 21.1299992 0.00000000 20.6280003 0.00000000 21.3969994 0.00000000 53544.9062 20.8169994 0.00000000 20.5650005 0.00000000 21.1459999 0.00000000
53544.9258 21.0049992 0.00000000 20.5020008 0.00000000 21.2719994 0.00000000 53544.9492 20.8799992 0.00000000 20.5020008 0.00000000 21.2089996 0.00000000 53544.9688 21.1930008 0.00000000 20.5020008 0.00000000 21.3339996 0.00000000
53544.9883 21.0680008 0.00000000 20.4389992 0.00000000 21.4589996 0.00000000 53545.0117 21.3810005 0.00000000 20.4389992 0.00000000 21.8980007 0.00000000 53545.0312 21.8190002 0.00000000 20.4389992 0.00000000 22.0230007 0.00000000
53545.0508 21.5680008 0.00000000 20.4389992 0.00000000 21.7730007 0.00000000 53545.0742 21.7560005 0.00000000 20.3759995 0.00000000 21.8349991 0.00000000 53545.0938 21.8810005 0.00000000 20.3759995 0.00000000 21.7730007 0.00000000
53545.1133 22.0070000 0.00000000 20.3759995 0.00000000 22.0230007 0.00000000 53545.1367 22.0070000 0.00000000 20.3759995 0.00000000 21.8980007 0.00000000 53545.1562 21.9440002 0.00000000 20.3759995 0.00000000 21.8980007 0.00000000
53545.1758 21.8810005 0.00000000 20.3759995 0.00000000 21.8980007 0.00000000 53545.1992 21.8810005 0.00000000 20.3759995 0.00000000 21.8349991 0.00000000 53545.2188 21.8810005 0.00000000 20.3759995 0.00000000 21.8349991 0.00000000
53545.2383 21.8810005 0.00000000 20.3759995 0.00000000 21.8349991 0.00000000 53545.2617 21.8190002 0.00000000 20.3759995 0.00000000 21.7730007 0.00000000 53545.2812 21.7560005 0.00000000 20.6910000 0.00000000 21.7730007 0.00000000
53545.3008 21.7560005 0.00000000 21.1310005 0.00000000 21.7730007 0.00000000 53545.3242 21.6940002 0.00000000 21.2570000 0.00000000 21.6469994 0.00000000 53545.3438 21.6940002 0.00000000 21.3199997 0.00000000 2
The matlab code used to generate “all_obs_ts.dat”
clc;clear;close all;
filename=('all_obs_ts.dat');
day=365;
start_date='01-Jan-2019 00:00:00';
end_date='31-Dec-2019 23:00:00';
obs_num=1;
pos_fvcom=6600;
% load mooring data
ncload('/Users/chuyanzhao/Desktop/gls_observation_data/obs_glerl_mooring/glerl_lake_superior_temperature_mooring_2018_2020.nc','sea_water_temperature','time');
temp=sea_water_temperature(:,[6:21]);
tmp_ind=temp<-1000;
temp(tmp_ind)=NaN;
[p,q]=size(temp);
% load glsea lst data
load('/Users/chuyanzhao/Desktop/fvcom_output/lake_superior/year_2019/code/processed_data/interpolated_obs_LowRes.mat','temp_obs');
for i=1:day
temp_obs_hour(24*(i-1)+1:i*24)=temp_obs(i,pos_fvcom);
end
temp_obs_lst=temp_obs_hour';
reference_date='1858-11-17 00:00:00';
tmp_time_obs=datetime(datestr(time/86400+datenum(1970,1,1)));
start_n=find(tmp_time_obs==start_date);
end_n=find(tmp_time_obs==end_date);
time_dat=datenum(tmp_time_obs-reference_date);
% define the first column of .dat
tmp_data(:,1)=time_dat(start_n:end_n,:);
[m,n]=size(tmp_data);
% define the second column of .dat (lst)
tmp_data(:,2)=temp_obs_lst(:);
% define the temperature of other layers
for i=4:2:(q+1)*2
tmp_data(:,i)=temp(start_n:end_n,i/2-1);
end
% define salnity
for i=3:2:(q+1)*2+1
tmp_data(:,i)=0.0;
end
% print
fileID = fopen(filename,'w');
output_format=strjoin([repmat({'%8.4f'},1,(q+1)*2+1),'\n']);
for i=1:obs_num
fprintf(fileID,'%d %d\n',i,m);
fprintf(fileID,output_format,tmp_data');
end
fclose(fileID);
- weight_control.dat
# the number of total mooring points
1
# mooring number; the fvcom grids number influenced by the mooring point
1 360
# node number; weight
5339 0.0197
5340 0.0276
5341 0.0195
5342 0.0140
5413 0.0238
5414 0.0228
5415 0.0259
5416 0.0356
5417 0.0526
5418 0.0354
5419 0.0258
5420 0.0187
5421 0.0129
5499 0.0365
5500 0.0388
5501 0.0436
5502 0.0568
5503 0.0882
5504 0.0915
5505 0.0577
5506 0.0445
5507 0.0336
5508 0.0250
5509 0.0154
5510 0.0109
5588 0.0339
5589 0.0522
5590 0.0589
5591 0.0678
5592 0.0793
5593 0.1297
5594 0.1395
5595 0.1443
5596 0.1045
5597 0.0719
5598 0.0569
5599 0.0449
5600 0.0377
5601 0.0245
The matlab code used to generate “weight_control.dat"
%% load data
clc;clear;close all;
% load fvcon grids
ncload(['/Users/chuyanzhao/Desktop/fvcom_output/lake_superior/year_2019/raw_data/da/low_res/sst/superior_0001.nc'],'lon','lat');
[m,n]=size(lat);
lat_mod=lat;
lon_mod=lon-360;
% load boundary
map=load('/Users/chuyanzhao/Desktop/FVCOM_output/Lake_Superior/boundary_superior.dat');
% load mooring location
ncload('/Users/chuyanzhao/Desktop/gls_observation_data/obs_glerl_mooring/glerl_lake_superior_temperature_mooring_2018_2020.nc','lat','lon');
lat_obs=lat;
lon_obs=lon;
obs_num=1;
% load glsea lst and fvcom lst to calculate error correlation
load('/Users/chuyanzhao/Desktop/fvcom_output/lake_superior/year_2019/code/processed_data/interpolated_obs_LowRes.mat','temp_mod','temp_obs');
%% define weight by error correlation
corr_bar=0.87;
dis_bar=84.0;
% find the nearest node of the mooring location
xy_node=[lat_mod lon_mod];
xy_obs=[lat_obs lon_obs; lat_obs lon_obs];
loc=find_nearest_group(xy_node,xy_obs);
pos=loc(1,1);
lst_mod=mean_by_segments(temp_mod,[0:24:8760],1);
% calculate the error correlation between every node and mooring location
for i=1:m
tmp_corr=corrcoef(lst_mod(:,pos)-temp_obs(:,pos),lst_mod(:,i)-temp_obs(:,i));
correlation(i)=tmp_corr(1,2);
end
% define weight
for i=1:m
distance2obs(i)=sqrt((lon_mod(i)-lon_obs)^2+(lat_mod(i)-lat_obs)^2)/180*pi*6371;
norm_tmp(i)=normpdf(distance2obs(i),0,dis_bar/3);
end
tmp_max=max(norm_tmp);
for i=1:m
if ((distance2obs(i)<dis_bar)&(correlation(i)>corr_bar))
weight_corr(i)=normpdf(distance2obs(i),0,dis_bar/3)/tmp_max*correlation(i);
else
weight_corr(i)=0.0;
end
end
% save no zero weight to data
k=1;
for i=1:m
if weight_corr(i)>0
data(k,1)=i;
data(k,2)=weight_corr(i);
k=k+1;
end
end
[size1_data,size2_data]=size(data);
%% print
filename1=['/Users/chuyanzhao/Desktop/fvcom_output/lake_superior/year_2019/code/processed_data/weight_control.dat'];
fileID1 = fopen(filename1,'w');
fprintf(fileID1,'%d\n',obs_num);
fprintf(fileID1,'%d %d\n',obs_num,size1_data);
fprintf(fileID1,'%d %8.4f\n',data');
fclose(fileID1);
%% figure the weight
figure(222);
clf;hold on;
cmin=0;cinc=0.1;cmax=1;
ncc=round([cmax-cmin]/cinc);
ccmap=flipud(spectral(ncc));
hh801=subplot(1,1,1);cla;hold on;box on;
scattercontourf(lon_mod,lat_mod,weight_corr',[cmin:cinc:cmax]);
scatter(lon_mod(pos),lat_mod(pos),'SizeData',150,'MarkerFaceColor','black');
text(lon_mod(pos),lat_mod(pos),'obs loc');
plot_fvcom_obc(map,[210 180 140]/255);
caxis([cmin cmax]);
colormap(hh801,ccmap);
lcb=colorbar('location','eastoutside');
扫一扫
73
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
