Commit 02754096 authored by Dany Dumont's avatar Dany Dumont
Browse files

README.md + formattage des scripts

parent f1cc58c8
SDM stands for Surface Drift Models that have been developped to improve the accuracy of floating object surface drift trajectory forecasts given wind and ocean currents. The package contains the code for four different drift models labelled A, B, C and D and that are described below.
#### 1. Control model
Model A is the most commonly used drift model where the surface drift is equal to the ocean current velocity to which is added a fraction of the 10-m wind velocity.
#### 2. Extrapolating surface currents
#### 5. References
Tamtare, T., D. Dumont, C. Chavanne (2019) Extrapolating eulerian currents for improving surface drift forecast, J. Oper. Oceanogr., submitted 28 September 2018, re-submitted 7 March 2019.
% MAIN_RK - Programme principal qui transporte des particules
% Inputs: Wind, current
%
% Outputs: particule position in lon and lat
% Other m-files required: RK.m, rotate_vector.m
% Author: michel tamtare
% ISMER-UQAR/polr
function [data] = main_RKmA(dat1,dat2,lat1,lon1)
% MAIN_RKMA - Programme principal qui transporte des particules ...
%
% Syntax: [data] = main_RKmA(dat1,dat2,lat1,lon1)
%
% Inputs:
% dat1 - date de debut de derive du spot
% dat2 - date de fin de derive du spot
% lat1 -
% lon1 -
%
% Outputs:
% data - particule position in lon and lat
%
% Example:
% [data] = main_RKmA(dat1,dat2,lat1,lon1)
%
% Other m-files required: RK.m, rotate_vector.m, ll2xpyp.m
% Subfunctions: none
% MAT-files required:
%
% See also: MAIN_RKMB, MAIN_RKMC, MAIN_RKMD
% Author: Michel Tamtare
% ISMER-UQAR
% email: tamkpanka.tamtare@uqar.ca
% juin 2018; Last revision:
% Last revision: June 2018
%------------- BEGIN CODE ------------------------------------------------------------------------------------------------------------------------------------
%*************************************************************************************************************************************************************
% Definition des constantes:
% Definition des constantes:
%*************************************************************************************************************************************************************
% h............................Pas de temps pour le transport de particules (1 = 1 tstep
......@@ -37,7 +49,6 @@
% k............................Indice de latitude dans la grille du modèle
% iii..........................Indice de longitude dans la grille du modèle
%*************************************************************************************************************************************************************
function [data] = main_RKmA(dat1,dat2,lat1,lon1)
close all
%clear mex
......@@ -45,205 +56,172 @@ close all
% Configuration
%*************************************************************************************************************************************************************
h = 3600; % Pas de temps (sec)
k = [];
iii = [];
reduit = 0; % Pour n'utiliser qu'une portion de la grille (plus rapide)
sens = 1; % -1 = backward
theta = -49; % Angle de rotation de la grille du golfe
dx = 5000; % Resolution de la grille (m)
if reduit ==1 % Modele reduit
%*********************** Cellules moullées *******************************************************************************************************************
i1 = 80;
i2 = 150;
k1 = 80;
k2 = 236;
h = 3600; % Pas de temps (sec)
k = [];
iii = [];
reduit = 0; % Pour n'utiliser qu'une portion de la grille (plus rapide)
sens = 1; % -1 = backward
theta = -49; % Angle de rotation de la grille du golfe
dx = 5000; % Resolution de la grille (m)
if reduit == 1 % Modele reduit
%*********************** Cellules mouillées ******************************************************************************************************************
i1 = 80;
i2 = 150;
k1 = 80;
k2 = 236;
else % Ensemble des cellules
i1 = 2;
i2 = 150;
k1 = 2;
k2 = 236;
i1 = 2;
i2 = 150;
k1 = 2;
k2 = 236;
end
[iii,k,status] = ll2xpyp(lon1,lat1,'gsl5km')
% trouver les indices des cellules mouill??es dans le domaine choisi
% trouver les indices des cellules mouillees dans le domaine choisi
t = 1;
for ii = i1:i2
for kk = k1:k2
%if ( ~isnan(test(ii,kk)) )
kcomp(t) = kk;
icomp(t) = ii;
t = t+1;
%end
end
for kk = k1:k2
%if ( ~isnan(test(ii,kk)) )
kcomp(t) = kk;
icomp(t) = ii;
t = t+1;
%end
end
end
%******************************************************************Fin de la configuration*********************************************************************
% ******************* Fin de la configuration *****************************
%****************************************************************Chargement des fichiers **********************************************************************
% ******************* Chargement des fichiers *****************************
%addpath('/share/archives/partage_lasso/ww3/gsl038_2012-2015_stokes/2013/gridded_dat')
%addpath('/home/tamt0001/these/data/drifter/2014/gridded_dat') ;
%% ************ Chargement des fichiers ***********************************************************************************************************************
addpath('/home/tamt0001/these/data/drifter/BAS/mat/GSL/datagsl') ;
% addpath('/home/tamt0001/these/data/2015/G2L') ;
% load 'G2L/G2Locean_2015.mat';
%addpath('/home/tamt0001/these/data/2015/G2L') ;
%load 'G2L/G2Locean_2015.mat';
load 'juin12014.mat';
dat= datO;
uc=u;
vc=v;
dat = datO;
uc = u;
vc = v;
%% vent du model GEM (test avec le spot 73)
load 'GEM2014.mat';
uw=Uw;
vw=Vw;
datw=datg;
uw = Uw;
vw = Vw;
datw = datg;
uw = permute(uw,[3 1 2]);
vw = permute(vw,[3 1 2]);
T1 = find(dat<=dat1+0.101); % % date1 est la date de d??but de d??rive du spot
T2 = find(dat<=dat2+0.001); % %date2 est la date de fin de d??rive du spot (voir tr_compare)
T1 = find(dat <= dat1 + 0.101); % date1 est la date de debut de derive du spot
T2 = find(dat <= dat2 + 0.001); % date2 est la date de fin de derive du spot (voir tr_compare)
T1=T1(end);
T1 = T1(end);
T2 = T2(end);
T2=T2(end);
T3 = find(datw <= dat1 + 0.001); % date1 est la date de debut de derive du spot
T4 = find(datw < dat2 + 0.001); % date2 est la date de fin de derive du spot (voir tr_compare)
T3 = find(datw<=dat1+0.001); % % date1 est la date de d??but de d??rive du spot
T4 = find(datw<dat2+0.001); % %date2 est la date de fin de d??rive du spot (voir tr_compare)
T3=T3(end);
T4=T4(end);
T3 = T3(end);
T4 = T4(end);
% keyboard
%**********************************************************************************************************************************************************
% 2. Transporter les particules avec Runge-Kutta 4
%**********************************************************************************************************************************************************
tsteps1 = T1:T2-1;
tsteps11=tsteps1+1;
tsteps2 = T3:T4-1;
tsteps1 = T1:T2-1;
tsteps11 = tsteps1+1;
tsteps22=tsteps2+1;
tsteps2 = T3:T4-1;
tsteps22 = tsteps2+1;
% keyboard
% tic
for tsteps = 1:length(tsteps1);
tsteps
% toc
% disp([ num2str(tsteps1-T1+1) '/' num2str(T2-T1) ])
uc1 = squeeze(uc(tsteps1(tsteps),:,:));
vc1 = squeeze(vc(tsteps1(tsteps),:,:));
uc2 = squeeze(uc(tsteps11(tsteps),:,:));
vc2 = squeeze(vc(tsteps11(tsteps),:,:));
uw1 = squeeze(uw(tsteps2((tsteps)),:,:));
vw1 = squeeze(vw(tsteps2(tsteps),:,:));
uw2 = squeeze(uw(tsteps22(tsteps),:,:));
vw2 = squeeze(vw(tsteps22(tsteps),:,:));
% keyboard
% ici la valeur de la fraction du vent en nombre complexe (modèleB)
mA= 0.0234 + sqrt(-1).*0.0015;
U1=uc1+sqrt(-1).*vc1;
Uw1=uw1+sqrt(-1).*vw1;
U2=uc2+sqrt(-1).*vc2;
Uw2=uw2+ sqrt(-1).*vw2;
Uc1=U1+mA.*Uw1;
Uc2=U2+mA.*Uw2;
% % keyboard
uc1=real(Uc1);
vc1=imag(Uc1);
uc2=real(Uc2);
vc2=imag(Uc2);
clear uw1 vw1 uw2 vw2
for t = 1:length(kcomp)
u1l(t) = squeeze(uc1(icomp(t),kcomp(t)));
v1l(t) = squeeze(vc1(icomp(t),kcomp(t)));
u2l(t) = squeeze(uc2(icomp(t),kcomp(t)));
v2l(t) = squeeze(vc2(icomp(t),kcomp(t)));
end
%keyboard
clear u1 v1 u2 v2
[u1,v1]=rotate_vector(u1l,v1l,theta);
[u2,v2]=rotate_vector(u2l,v2l,theta);
% keyboard
clear u1l v1l u2l v2l
% u1(isnan(u1)) = 0;
% v1(isnan(v1)) = 0;
% u2(isnan(u2)) = 0;
% v2(isnan(v2)) = 0;
%*******************************************************************************************************************************************************
% Fonction RK.m pour le transport
%******************************************************************************************************************************************************
umid = ( u1 + u2 )/2;
vmid = ( v1 + v2 )/2;
[k,iii] = RK(u1,u2,-v1,-v2,umid,-vmid,k,iii,h,dx,kcomp,icomp);
% toc
% disp([ num2str(tsteps1-T1+1) '/' num2str(T2-T1) ])
uc1 = squeeze(uc(tsteps1 (tsteps),:,:));
vc1 = squeeze(vc(tsteps1 (tsteps),:,:));
uc2 = squeeze(uc(tsteps11(tsteps),:,:));
vc2 = squeeze(vc(tsteps11(tsteps),:,:));
uw1 = squeeze(uw(tsteps2 (tsteps),:,:));
vw1 = squeeze(vw(tsteps2 (tsteps),:,:));
uw2 = squeeze(uw(tsteps22(tsteps),:,:));
vw2 = squeeze(vw(tsteps22(tsteps),:,:));
% keyboard
% ici la valeur de la fraction du vent en nombre complexe (modèleB)
mA = 0.0234 + sqrt(-1).*0.0015;
U1 = uc1+sqrt(-1).*vc1;
Uw1 = uw1+sqrt(-1).*vw1;
U2 = uc2+sqrt(-1).*vc2;
Uw2 = uw2+ sqrt(-1).*vw2;
Uc1 = U1+mA.*Uw1;
Uc2 = U2+mA.*Uw2;
% % keyboard
uc1 = real(Uc1);
vc1 = imag(Uc1);
uc2 = real(Uc2);
vc2 = imag(Uc2);
clear uw1 vw1 uw2 vw2
for t = 1:length(kcomp)
u1l(t) = squeeze(uc1(icomp(t),kcomp(t)));
v1l(t) = squeeze(vc1(icomp(t),kcomp(t)));
u2l(t) = squeeze(uc2(icomp(t),kcomp(t)));
v2l(t) = squeeze(vc2(icomp(t),kcomp(t)));
end
%keyboard
clear u1 v1 u2 v2
[u1,v1] = rotate_vector(u1l,v1l,theta);
[u2,v2] = rotate_vector(u2l,v2l,theta);
% keyboard
clear u1l v1l u2l v2l
% u1(isnan(u1)) = 0;
% v1(isnan(v1)) = 0;
% u2(isnan(u2)) = 0;
% v2(isnan(v2)) = 0;
%*******************************************************************************************************************************************************
% Fonction RK.m pour le transport
%******************************************************************************************************************************************************
umid = ( u1 + u2 )/2;
vmid = ( v1 + v2 )/2;
[k,iii] = RK(u1,u2,-v1,-v2,umid,-vmid,k,iii,h,dx,kcomp,icomp);
end
%time = (dat1:1/24:dat2);
time = dat(T1:T2);
[lon,lat,status] = xpyp2ll(iii,k,'gsl5km');
%keyboard
%save(sprintf('path/%s.mat',name),'time','iii','k','lon','lat');
if numel(time) ~= numel(lat)
disp('numel(time) ~= numel(lat)')
whos spot time lat lon
data = [NaN];
return
disp('numel(time) ~= numel(lat)')
whos spot time lat lon
data = [NaN];
return
else
data = [time,lon,lat];
data = [time,lon,lat];
end
% MAIN_RK - Programme principal qui transporte des particules
% Inputs: Wind, current
%
% Outputs: particule position in lon and lat
% Other m-files required: RK.m, rotate_vector.m
% Author: michel tamtare
% ISMER-UQAR/polr
function [data] = main_RKmB(dat1,dat2,lat1,lon1)
% MAIN_RKMB - Programme principal qui transporte des particules ...
%
% Syntax: [data] = main_RKmB(dat1,dat2,lat1,lon1)
%
% Inputs:
% dat1 - date de debut de derive du spot
% dat2 - date de fin de derive du spot
% lat1 -
% lon1 -
%
% Outputs:
% data - particule position in lon and lat
%
% Example:
% [data] = main_RKmB(dat1,dat2,lat1,lon1)
%
% Other m-files required: RK.m, rotate_vector.m, ll2xpyp.m
% Subfunctions: none
% MAT-files required:
%
% See also: MAIN_RKMA, MAIN_RKMC, MAIN_RKMD
% Author: Michel Tamtare
% ISMER-UQAR
% email: tamkpanka.tamtare@uqar.ca
% juin 2018; Last revision:
% Last revision: June 2018
%------------- BEGIN CODE ------------------------------------------------------------------------------------------------------------------------------------
%*************************************************************************************************************************************************************
% Definition des constantes:
% Definition des constantes:
%*************************************************************************************************************************************************************
% h............................Pas de temps pour le transport de particules (1 = 1 tstep
......@@ -37,7 +49,7 @@
% k............................Indice de latitude dans la grille du modèle
% iii..........................Indice de longitude dans la grille du modèle
%*************************************************************************************************************************************************************
function [data] = main_RKmB(dat1,dat2,lat1,lon1)
close all
%clear mex
......@@ -45,25 +57,25 @@ close all
% Configuration
%*************************************************************************************************************************************************************
h = 3600; % Pas de temps (sec)
k = [];
iii = [];
h = 3600; % Pas de temps (sec)
k = [];
iii = [];
reduit = 0; % Pour n'utiliser qu'une portion de la grille (plus rapide)
sens = 1; % -1 = backward
theta = -49; % Angle de rotation de la grille du golfe
dx = 5000; % Resolution de la grille (m)
if reduit ==1 % Modele reduit
%*********************** Cellules moullées *******************************************************************************************************************
i1 = 80;
i2 = 150;
k1 = 80;
k2 = 236;
sens = 1; % -1 = backward
theta = -49; % Angle de rotation de la grille du golfe
dx = 5000; % Resolution de la grille (m)
if reduit == 1 % Modele reduit
%*********************** Cellules moullées *******************************************************************************************************************
i1 = 80;
i2 = 150;
k1 = 80;
k2 = 236;
else % Ensemble des cellules
i1 = 2;
i2 = 150;
k1 = 2;
k2 = 236;
i1 = 2;
i2 = 150;
k1 = 2;
k2 = 236;
end
......@@ -72,17 +84,17 @@ end
% trouver les indices des cellules mouill??es dans le domaine choisi
t = 1;
for ii = i1:i2
for kk = k1:k2
%if ( ~isnan(test(ii,kk)) )
kcomp(t) = kk;
icomp(t) = ii;
t = t+1;
%end
end
for kk = k1:k2
%if ( ~isnan(test(ii,kk)) )
kcomp(t) = kk;
icomp(t) = ii;
t = t+1;
%end
end
end
......@@ -111,26 +123,26 @@ datw=datg;
uw = permute(uw,[3 1 2]);
vw = permute(vw,[3 1 2]);
%% Ici nous faisons une extrapolation du courant en surface pour en faire une prévision de dérive en surface
%% Extrapollation linéaire
%% Extrapollation linéaire
z1= 2.5 ; z2=7.5 ;
z1= 2.5 ; z2=7.5 ;
uc=1./2*(uca+ucb)-1./2*((ucb-uca)./(z2-z1))*(z1+z2);
vc=1./2*(vca+vcb)-1./2*((vcb-vca)./(z2-z1))*(z1+z2);
% keyboard
T1 = find(dat<=dat1+0.101); % % date1 est la date de d??but de d??rive du spot
T2 = find(dat<=dat2+0.001); % %date2 est la date de fin de d??rive du spot (voir tr_compare)
T1 = find(dat<=dat1+0.101); % % date1 est la date de d??but de d??rive du spot
T2 = find(dat<=dat2+0.001); % %date2 est la date de fin de d??rive du spot (voir tr_compare)
T1=T1(end);
T2=T2(end);
T1=T1(end);
T2=T2(end);
T3 = find(datw<=dat1+0.001); % % date1 est la date de d??but de d??rive du spot
T4 = find(datw<dat2+0.001); % %date2 est la date de fin de d??rive du spot (voir tr_compare)
T3=T3(end);
T4=T4(end);
T3 = find(datw<=dat1+0.001); % % date1 est la date de d??but de d??rive du spot
T4 = find(datw<dat2+0.001); % %date2 est la date de fin de d??rive du spot (voir tr_compare)
T3=T3(end);
T4=T4(end);
% keyboard
%****************************************************************************************************************************************************
......@@ -148,89 +160,78 @@ tsteps22=tsteps2+1;
% tic
for tsteps = 1:length(tsteps1);
tsteps
% toc
% disp([ num2str(tsteps1-T1+1) '/' num2str(T2-T1) ])
uc1 = squeeze(uc(tsteps1(tsteps),:,:));
vc1 = squeeze(vc(tsteps1(tsteps),:,:));
uc2 = squeeze(uc(tsteps11(tsteps),:,:));
vc2 = squeeze(vc(tsteps11(tsteps),:,:));
uw1 = squeeze(uw(tsteps2((tsteps)),:,:));
vw1 = squeeze(vw(tsteps2(tsteps),:,:));
uw2 = squeeze(uw(tsteps22(tsteps),:,:));
vw2 = squeeze(vw(tsteps22(tsteps),:,:));
% keyboard
% ici la valeur de la fraction du vent en nombre complexe (modèleB)
mB= 0.0150 + sqrt(-1).*0.0005;
U1=uc1+sqrt(-1).*vc1;
Uw1=uw1+sqrt(-1).*vw1;
U2=uc2+sqrt(-1).*vc2;
Uw2=uw2+ sqrt(-1).*vw2;
Uc1=U1+mB.*Uw1;
Uc2=U2+mB.*Uw2;
% % keyboard
uc1=real(Uc1);
vc1=imag(Uc1);
uc2=real(Uc2);
vc2=imag(Uc2);
% % save(sprintf('drifter_mat2/U.mat'),'datO','uc12','vc12','uc22','vc22')
% % keyboard
clear uw1 vw1 uw2 vw2
for t = 1:length(kcomp)
u1l(t) = squeeze(uc1(icomp(t),kcomp(t)));
v1l(t) = squeeze(vc1(icomp(t),kcomp(t)));
u2l(t) = squeeze(uc2(icomp(t),kcomp(t)));
v2l(t) = squeeze(vc2(icomp(t),kcomp(t)));
end
%keyboard
clear u1 v1 u2 v2
[u1,v1]=rotate_vector(u1l,v1l,theta);
[u2,v2]=rotate_vector(u2l,v2l,theta);
% keyboard
clear u1l v1l u2l v2l
% u1(isnan(u1)) = 0;
% v1(isnan(v1)) = 0;
% u2(isnan(u2)) = 0;
% v2(isnan(v2)) = 0;
umid = ( u1 + u2 )/2;
vmid = ( v1 + v2 )/2;
% toc
% disp([ num2str(tsteps1-T1+1) '/' num2str(T2-T1) ])
uc1 = squeeze(uc(tsteps1 (tsteps),:,:));