-
Notifications
You must be signed in to change notification settings - Fork 3
/
hycom2fvcom_iniTS_create.m
113 lines (100 loc) · 3.22 KB
/
hycom2fvcom_iniTS_create.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
%==========================================================================
% HYCOM2FVCOM:
% Create TS initial file from HYCOM datset
%
% input :
% fgrid --- the file containing the FVCOM grid (nc, grd, or 2dm)
% fhycom --- hycom data (NetCDF)
% fout --- output initial TS file name and path (NetCDF)
%
% output :
% initial TS file
%
% Siqi Li, SMAST
% 2022-11-30
%
% Updates:
% 2023-01-03 Siqi Li Convert the fvcom longitude into [0 360]
% 2023-01-04 Siqi Li Used new interpolating method
%==========================================================================
% addpath('~/tools/matFVCOM')
% addpath('~/tools/matFigure')
clc
clear
%--------------------------------------------------------------------------
% Input
fgrid = '../data/Gom7_v5_nest_CascoBay5.1.2dm';
fhycom = '../hycom/hycom_20161221_0000.nc';
fout = '../output/gom7_ini_ts_20161221.nc';
%--------------------------------------------------------------------------
% Read FVCOM grid and sigma
f = f_load_grid(fgrid);
[lon_fvcom, lat_fvcom] = sp_proj('1802', 'inverse', f.x, f.y, 'm');
% Read HYCOM data
lon0 = ncread(fhycom, 'lon');
lat0 = ncread(fhycom, 'lat');
depth0 = ncread(fhycom, 'depth');
t0 = ncread(fhycom, 'water_temp');
s0 = ncread(fhycom, 'salinity');
time0 = ncread(fhycom, 'time')/24 + datenum(2000,1,1);
% Dimensions
nz0 = length(depth0);
% Interpolation
t = nan(f.node, nz0);
s = nan(f.node, nz0);
disp('Calculating interpolating weight.')
wh = interp_2d_calc_weight('GLOBAL_BI', lon0, lat0, lon_fvcom, lat_fvcom);
for iz = 1 : nz0
disp(['Interpolating the ' num2str(iz) 'th layer of ' num2str(nz0) ' layers.'])
% Horizontal interpolation
t_layer = interp_2d_via_weight(t0(:,:,iz), wh);
s_layer = interp_2d_via_weight(s0(:,:,iz), wh);
% Fill the points that are recognize as land in HYCOM
t_layer = f_fill_missing(f, t_layer);
s_layer = f_fill_missing(f, s_layer);
% Set the layer below the depth as nan
if iz> 2
k_bot = f.h <= depth0(iz);
t_layer(k_bot) = nan;
s_layer(k_bot) = nan;
end
% Store the data
t(:,iz) = t_layer;
s(:,iz) = s_layer;
% % t
% t_layer = t0(:,:,iz);
% kt = ~isnan(t_layer);
% Ft = scatteredInterpolant(xx0(kt), yy0(kt), t_layer(kt), 'linear', 'nearest');
% t(:,iz) = Ft(x360, f.y);
% % s
% s_layer = s0(:,:,iz);
% ks = ~isnan(s_layer);
% Fs = scatteredInterpolant(xx0(kt), yy0(kt), s_layer(kt), 'linear', 'nearest');
% s(:,iz) = Fs(x360, f.y);
end
% Fill the points under the bathemetry
t = fillmissing(t, 'nearest', 2);
s = fillmissing(s, 'nearest', 2);
% Plot (Uncomment this part for figures)
figdir = '../output/';
close all
for iz = 1 : nz0
disp(['----' num2str(iz, '%2.2d')])
figure('Visible', 'off')
hold on
subplot(2,1,1)
f_2d_image(f, t(:,iz));
colorbar
caxis([0 35])
title([num2str(depth0(iz)) ' m (' num2str(iz, '%2.2d') '): Temperature'])
subplot(2,1,2)
f_2d_image(f, s(:,iz));
colorbar
caxis([25 40])
title([num2str(depth0(iz)) ' m (' num2str(iz, '%2.2d') '): Salinity'])
ffig = [figdir '/init_ts_' datestr(time0, 'yyyymmddHH') '_layer_' num2str(iz,'%2.2d') '.png'];
mf_save(ffig);
close
end
% Write initial TS output
write_initial_ts(fout, -depth0, t, s, time0);