-
Notifications
You must be signed in to change notification settings - Fork 0
/
ga_advection_ariane.m
131 lines (116 loc) · 6.41 KB
/
ga_advection_ariane.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
function positions=ga_advection_ariane(pos_ini,name_curr,varargin)
%% GA_ADVECTION_ARIANE: computes 2D surface current trajectories using Ariane
%
% Use:
% positions=ga_advection_ariane(pos_ini,name_curr,varargin)
%
% Required inputs:
% pos_ini initial positions, as a matrix with rows [lon,lat]
% name_curr prefix of the current dataset to be used (netcdf files as generated by utils/ga_write_ariane_currents into dir_ariane_global/currents_data)
%
% Optional inputs:
% 'dt' output resolution (default 0.2 days). NOTE: for Ariane dt impacts the output resolution only, NOT the result (small dt slightly increases processing time).
% 'time0' initial time (default first current time step)
% 'nbdays_advec' trajectory length (in days, default maximum ie as long as currents are available)
% 'backward' backward instead of forward calculations
%
% For instance, to change nbdays_advec:
% positions=ga_advection_ariane(pos_ini,name_curr,'nbdays_advec',100)
% To run a backward calculation starting at time0 (which is then the end of the trajectory):
% positions=ga_advection_ariane(pos_ini,name_curr,'time0',time0,'backward')
%
% Monique Messié, 2021 for public version
% Reference: Messié, M., Petrenko, A., Doglioli, A. M., Aldebert, C., Martinez, E., Koenig, G., Bonnet, S., and Moutin, T. (2020)
% The delayed island mass effect: How islands can remotely trigger blooms in the oligotrophic ocean.
% Geophysical Research Letters, 47(2), e2019GL085282, https://doi.org/10.1029/2019GL085282
global dir_ariane_global
[arg,flag]=ga_read_varargin(varargin,{'dt',0.2,'time0',[],'nbdays_advec',[]},{'backward'});
if isempty(pos_ini)
positions=struct();
for varname={'lon2D','lat2D','time','lon_ini','lat_ini'}, varname=varname{:}; positions.(varname)=NaN; return, end
end
% load the current grid and time
rep_curr=[dir_ariane_global,'currents_data/'];
ncfile=[rep_curr,'meshmask_',name_curr,'.nc'];
gridcurr=struct();
gridcurr.u=struct(); gridcurr.u.lon=unique(double(ncread(ncfile,'xu'))); gridcurr.u.lat=unique(double(ncread(ncfile,'yu')));
gridcurr.v=struct(); gridcurr.v.lon=unique(double(ncread(ncfile,'xv'))); gridcurr.v.lat=unique(double(ncread(ncfile,'yv')));
if min(gridcurr.u.lon(:))>max(pos_ini(:,1))
shift_lon_factor=360;
elseif max(gridcurr.u.lon(:))<min(pos_ini(:,1))
shift_lon_factor=-360;
else, shift_lon_factor=0;
end
list_curr=ga_dir2filenames(rep_curr,[name_curr,'_0*.nc']);
time_curr=ones(length(list_curr),1)*NaN;
for itime=1:length(list_curr), time_curr(itime)=ncread([rep_curr,list_curr{itime}],'time'); end
% write namelist
if isempty(arg.time0)
if flag.backward, arg.time0=max(time_curr); else, arg.time0=min(time_curr); end
end
if arg.time0<min(time_curr) || arg.time0>max(time_curr), error('time0 must be within the current time period'), end
if isempty(arg.nbdays_advec), arg.nbdays_advec=max(time_curr)-arg.time0; end
if flag.backward, itime_curr=[interp1(time_curr,1:length(time_curr),arg.time0-arg.nbdays_advec,'nearest') interp1(time_curr,1:length(time_curr),arg.time0,'nearest')];
else, itime_curr=[interp1(time_curr,1:length(time_curr),arg.time0,'nearest') interp1(time_curr,1:length(time_curr),arg.time0+arg.nbdays_advec,'nearest')];
end
if isnan(itime_curr(2)), disp('Not enough current temporal coverage, adjusting nbdays_advec')
itime_curr(2)=length(time_curr); nbdays_advec=time_curr(end)-arg.time0;
else, nbdays_advec=arg.nbdays_advec;
end
ga_write_ariane_namelist(name_curr,'dt',arg.dt,'nbdays_advec',nbdays_advec,...
'dtcurr',time_curr(2)-time_curr(1),'itime_curr',itime_curr,varargin{:});
% genere posiition matrix (x, y, z, t0, weight) using U-V grids
pos_ini_xy=pos_ini*NaN;
for ipts=1:size(pos_ini,1)
pos_ini_xy(ipts,:)=[interp1(gridcurr.u.lon,1:length(gridcurr.u.lon),pos_ini(ipts,1)+shift_lon_factor),...
interp1(gridcurr.v.lat,1:length(gridcurr.v.lat),pos_ini(ipts,2))];
end
iround=round(pos_ini_xy)==pos_ini_xy; pos_ini_xy(iround)=pos_ini_xy(iround)+1E-3; % slightly shifts some positions to avoid an issue with Ariane
it0=interp1(time_curr,1:length(time_curr),arg.time0,'nearest')-itime_curr(1)+1; % this is in Ariane format, ie indexed on the re-indexed current files (depends on dtcurr)
iok_initpositions=~max(isnan(pos_ini_xy),[],2);
mat_positions_file=[pos_ini_xy(iok_initpositions,:),repmat([-1 it0 1],sum(iok_initpositions),1)];
% write positions file
posfile=[dir_ariane_global,'initial_positions.txt'];
fid=fopen(posfile,'w');
for iline=1:size(mat_positions_file,1)
fprintf(fid,'%.3f %.3f %.3f %.3f %.3f\n',mat_positions_file(iline,:));
end
fclose(fid);
% run ariane - see https://www.mathworks.com/matlabcentral/answers/570193-netcdf-error-in-matlab2019a
oldPath = getenv('LD_LIBRARY_PATH');
cleanVar = onCleanup(@()setenv('LD_LIBRARY_PATH', oldPath));
setenv('LD_LIBRARY_PATH','/usr/lib/x86_64-linux-gnu/')
unix(['cd ',dir_ariane_global,'; ariane > ariane_output.txt;']);
clear('cleanVar')
% load results: time (and check it is what we expected)
time_expected=(0:arg.dt:nbdays_advec)';
ncfile=[dir_ariane_global,'ariane_trajectories_qualitative.nc'];
time_all=ncread(ncfile,'traj_time'); time_all=time_all(1,:)';
time_all(time_all==ncreadatt(ncfile,'traj_time','missing_value'))=NaN;
if flag.backward, time_all=-time_all; time_expected=-time_expected; end
if max(abs(time_all(:)-time_expected(:)))>1E-10, error('check time expected'), end
% load results: positions
positions=struct();
positions.time=time_expected+arg.time0; % using time_expected because time_all is at 0 for particles outside the domain which screws up the time sorting in backward
positions.lon_ini=pos_ini(:,1);
positions.lat_ini=pos_ini(:,2);
positions.lon2D=nan(length(positions.lon_ini),arg.nbdays_advec/arg.dt+1);
positions.lat2D=nan(length(positions.lon_ini),arg.nbdays_advec/arg.dt+1);
itime_ariane=1:nbdays_advec/arg.dt+1;
if length(positions.time)>1
positions.lon2D(iok_initpositions,itime_ariane)=ncread(ncfile,'traj_lon')-shift_lon_factor;
positions.lat2D(iok_initpositions,itime_ariane)=ncread(ncfile,'traj_lat');
else
positions.lon2D(:,itime_ariane)=ncread(ncfile,'traj_lon')-shift_lon_factor;
positions.lat2D(:,itime_ariane)=ncread(ncfile,'traj_lat');
end
positions.lon2D(positions.lon2D==ncreadatt(ncfile,'traj_lon','missing_value'))=NaN;
positions.lat2D(positions.lat2D==ncreadatt(ncfile,'traj_lat','missing_value'))=NaN;
% Put positions back in the correct order for backward
if flag.backward
[~,isort]=sort(positions.time);
positions.time=positions.time(isort);
positions.lon2D=positions.lon2D(:,isort);
positions.lat2D=positions.lat2D(:,isort);
end
return