FVCOM & DA

Alice

已于 2023-10-06 03:01:03 修改

阅读量152

点赞数

文章标签： matlab

于 2023-09-09 05:56:43 首次发布

本文链接：https://blog.youkuaiyun.com/weixin_42222933/article/details/132729064

版权

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');