forked from rafacsantana/scripts_ww3
-
Notifications
You must be signed in to change notification settings - Fork 0
/
write_ccam_wind_nc4ww3.m
227 lines (191 loc) · 8.41 KB
/
write_ccam_wind_nc4ww3.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
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
clear
close all
run('/scale_wlg_persistent/filesets/home/santanarc/scripts/niwa/matlab/startup.m')
% Daily time
time_lima=datenum(1985,1,1,0,0,0):1:datenum(1990,1,1,0,0,0); % Graham Harrington ECAN request on 18/09/2023
%time_lima=datenum(1990,1,1,0,0,0):1:datenum(2014,12,31,0,0,0); % Graham Harrington ECAN request on 18/09/2023
%time_lima=datenum(1992,2,29,0,0,0):1:datenum(1992,2,29,0,0,0); % Graham Harrington ECAN request on 18/09/2023
time_lima=datenum(2014,01,01,0,0,0):1:datenum(2014,12,31,0,0,0); % Graham Harrington ECAN request on 18/09/2023
%time_lima=datenum(2015,01,01,0,0,0):1:datenum(2015,12,31,0,0,0); % Graham Harrington ECAN request on 18/09/2023
time_lima=datenum(2029,12,1,0,0,0):1:datenum(2030,12,31,0,0,0); % Graham Harrington ECAN request on 18/09/2023
time_lima=datenum(2031,1,1,0,0,0):1:datenum(2100,12,30,0,0,0); % Graham Harrington ECAN request on 18/09/2023
%/scale_wlg_nobackup/filesets/nobackup/niwa03150/WAVE/projections/historical_proc/6hrly_global/wind/2014/wind_historical_GFDL-ESM4_CCAM_6hrly_Global_raw_2014121700.nc
%Index exceeds the number of array elements (0).
%
%Error in write_wind_nc4ww3 (line 60)
% ua=squeeze((ncread(fname,'ua', [1 1 1 it(1)],[Inf Inf Inf it(end)-it(1)+1])));
% first merge u and v files
% cdo merge ua_historical_GFDL-ESM4_CCAM_6hrly_Global_raw.nc va_historical_GFDL-ESM4_CCAM_6hrly_Global_raw.nc wind_historical_GFDL-ESM4_CCAM_6hrly_Global_raw.nc
path_ou='/scale_wlg_nobackup/filesets/nobackup/niwa03150/WAVE/projections/historical_proc/6hrly_global/';
vars=[2];
vnames ={'ua/','sic/'}; % (ua for wind) always with /
if time_lima(1)<datenum(2015,1,1)
prefixs ={'ua_historical_','sic_historical_'}; % prior to 2015/01/01
path_in='/scale_wlg_nobackup/filesets/nobackup/niwa03150/WAVE/projections/historical/6hrly_global/';
else
prefixs ={'ua_ssp370_','sic_ssp370_'};
path_in='/scale_wlg_nobackup/filesets/nobackup/niwa03150/WAVE/projections/ssp370/6hrly_global/';
end
midfixs ={'GFDL-ESM4_CCAM_' ,'GFDL-ESM4_CCAM_'};
sufixs ={'6hrly_Global_raw','6hrly_Global_raw'};
scrsz=[1 1 1366 768];
scrsz=[2 42 958 953];
%scrsz=get(0,'screensize');
for vr=vars
vname ={vnames{vr}};
prefix={prefixs{vr}};
midfix={midfixs{vr}};
sufix ={sufixs{vr}};
fname=[path_in,vname{1},prefix{1},midfix{1},sufix{1},'.nc'];
display(['Loading: ',fname]);
% Loop in model time
for t=time_lima
vn=vname{1};
ptime=datestr(t,'YYYY/');
ftime=datestr(t,'YYYYmmDDHH');
if t==time_lima(1)
display(['Reading lon, lat, time']); tic;
FillValue=ncreadatt(fname,vn(1:end-1),'_FillValue');
lon=(ncread(fname,'lon'));
lat=(ncread(fname,'lat'));
timew=squeeze(double(ncread(fname,'time')));
if strcmp(prefix{1},'ua_historical_') || strcmp(prefix{1},'sic_historical_')
time_ini=datenum(1959,1,1);
else
time_ini=datenum(2015,1,1);
end
time_mat=timew./(24*60)+time_ini;
timed=datevec(time_mat);
toc
% climate models dont have leap years. removing 29/2 and adding a day to get to 31/21 of 2015/2100
% ww3 will read 29/02 as 28/02
timeww=timew;
time_matt=timeww./(24*60)+time_ini;
timedd=datevec(time_matt);
for i=1:length(timeww)
tt=time_matt(i);
td=datevec(tt);
if td(:,2)==2 & td(:,3)==29
it=find(timedd(:,1)==td(:,1) & timedd(:,2)==td(:,2) & timedd(:,3)==td(:,3));
timeww(it:end)=timeww(it:end)+24*60; % plus a day (minutes)
time_matt=timeww./(24*60)+time_ini;
timedd=datevec(time_matt);
end
end
time_matt=timeww./(24*60)+time_ini;
timedd=datevec(time_matt);
timed=timedd;
time_mat=time_matt;
timew=timeww;
end
% find indexes of selected day on nc file
td=datevec(t);
%it=find(timed(:,1)==td(:,1) & timed(:,2)==td(:,2) & timed(:,3)==td(:,3));
[dif it]=nanmin(abs(time_mat-t));
if strcmp(sufix{1},'6hrly_Global_raw')
it=it:it+4;
else % hourly
it=it:it+24;
end
% time that will be written
timewit=(t:1/(length(it)-1):t+1)-time_ini;
timewit=timewit(1:end)*(24*60);
%display(['Processing times: ']) % ,num2str(timew(it)./(24*60)+time_ini)])
%datevec(timewit./(24*60)+time_ini)
% need to create data for 2015/01/01 00:00
%if td(:,1)==2014 & td(:,2)==12 & td(:,3)==31
%end
if strcmp(prefix{1},'ua_historical_') || strcmp(prefix{1},'ua_ssp370_')
display(['Reading winds']); tic;
fname=[path_in,vname{1},prefix{1},midfix{1},sufix{1},'.nc'];
%fname=[path_in,'ua/','ua_historical_',midfix{1},sufix{1},'.nc'];
ua=squeeze((ncread(fname,'ua', [1 1 1 it(1)],[Inf Inf Inf it(end)-it(1)+1])));
if strcmp(prefix{1},'ua_historical_')
fname=[path_in,'va/','va_historical_',midfix{1},sufix{1},'.nc'];
elseif strcmp(prefix{1},'ua_ssp370_')
fname=[path_in,'va/','va_ssp370_',midfix{1},sufix{1},'.nc'];
end
va=squeeze((ncread(fname,'va', [1 1 1 it(1)],[Inf Inf Inf it(end)-it(1)+1])));
path_out=[path_ou,'wind/',ptime];
system(['mkdir -p ',path_out]);
fda=[path_out,'wind_historical_',midfix{1},sufix{1},'_',ftime,'.nc'];
display(['Saving: ',fda]);
if exist(fda)==2
system(['rm -rf ',fda]);
end
% creating nc and variables
ncid = netcdf.create(fda,'NETCDF4');
% time
dtime = netcdf.defDim(ncid,'time',length(it));
varid = netcdf.defVar(ncid,'time','double',dtime);
netcdf.putVar(ncid,varid,timewit)
netcdf.putAtt(ncid,varid,'standard_name',"time");
netcdf.putAtt(ncid,varid,'calendar',"gregorian");
if strcmp(prefix{1},'ua_historical_')
netcdf.putAtt(ncid,varid,'units',"minutes since 1959-01-01 00:00:00");
elseif strcmp(prefix{1},'ua_ssp370_')
netcdf.putAtt(ncid,varid,'units',"minutes since 2015-01-01 00:00:00");
end
% lon
dlon = netcdf.defDim(ncid,'lon',length(lon));
varid = netcdf.defVar(ncid,'lon','double',dlon);
netcdf.putVar(ncid,varid,lon)
% lat
dlat = netcdf.defDim(ncid,'lat',length(lat));
varid = netcdf.defVar(ncid,'lat','double',dlat);
netcdf.putVar(ncid,varid,lat)
% ua
uvar = netcdf.defVar(ncid,'ua','float', [dlon dlat dtime]);
netcdf.defVarFill(ncid,uvar,false,FillValue);
netcdf.putVar(ncid,uvar,ua)
netcdf.putAtt(ncid,uvar,'standard_name',"eastward_wind");
netcdf.putAtt(ncid,uvar,'units',"m s-1");
% va
vvar = netcdf.defVar(ncid,'va','float', [dlon dlat dtime]);
netcdf.defVarFill(ncid,vvar,false,FillValue);
netcdf.putVar(ncid,vvar,va)
netcdf.putAtt(ncid,vvar,'standard_name',"northward_wind");
netcdf.putAtt(ncid,vvar,'units',"m s-1");
netcdf.close(ncid);
elseif strcmp(prefix{1},'sic_historical_') || strcmp(prefix{1},'sic_ssp370_')
display(['Reading sic']); tic;
sic=squeeze((ncread(fname,'sic', [1 1 it(1)],[Inf Inf it(end)-it(1)+1])));
sic=sic./100;
path_out=[path_ou,vname{1},ptime];
system(['mkdir -p ',path_out]);
fda=[path_out,prefix{1},midfix{1},sufix{1},'_',ftime,'.nc'];
display(['Saving: ',fda]);
if exist(fda)==2
system(['rm -rf ',fda]);
end
% creating nc and variables
ncid = netcdf.create(fda,'NETCDF4');
% time
dtime = netcdf.defDim(ncid,'time',length(it));
varid = netcdf.defVar(ncid,'time','double',dtime);
netcdf.putVar(ncid,varid,timewit)
netcdf.putAtt(ncid,varid,'standard_name',"time");
netcdf.putAtt(ncid,varid,'calendar',"gregorian");
if strcmp(prefix{1},'sic_historical_')
netcdf.putAtt(ncid,varid,'units',"minutes since 1959-01-01 00:00:00");
elseif strcmp(prefix{1},'sic_ssp370_')
netcdf.putAtt(ncid,varid,'units',"minutes since 2015-01-01 00:00:00");
end
% lon
dlon = netcdf.defDim(ncid,'lon',length(lon));
varid = netcdf.defVar(ncid,'lon','double',dlon);
netcdf.putVar(ncid,varid,lon)
% lat
dlat = netcdf.defDim(ncid,'lat',length(lat));
varid = netcdf.defVar(ncid,'lat','double',dlat);
netcdf.putVar(ncid,varid,lat)
% sic
sicvar = netcdf.defVar(ncid,'sic','float', [dlon dlat dtime]);
netcdf.defVarFill(ncid,sicvar,false,FillValue);
netcdf.putVar(ncid,sicvar,sic)
netcdf.putAtt(ncid,sicvar,'standard_name',"sea_ice_area_fraction");
netcdf.putAtt(ncid,sicvar,'units',"1");
netcdf.close(ncid);
end % if strcmp(prefix
end
end