Skip to content
GitLab
Commit f1cc58c8 authored by Dany Dumont's avatar Dany Dumont
Browse files

Premier depot

parents
Hide whitespace changes
Inline Side-by-side
geos.m 0 → 100644
View file @ f1cc58c8
function [datO,LON,LAT,ugeo,vgeo]=geos(datO,wlev,lon,lat)
%% calcul du courant geostrophique
dat= datO;
ssh=wlev;
LON=lon;
LAT=lat;
g=9.81;
f=10^-4;
Psi=(ssh.*g)./f;
Psi = permute(Psi,[2 3 1]);
[m n p]=size(Psi);
%DT is in hour!!!
%% Elimination des marées
for i=1:m
for j=1:n
int=1; % 1h en heur
I = find(~isnan(Psi(i,j,:)));
if ~isempty(I)
[NAME,FREQ,TIDECON,XOUT(i,j,:)]=t_tide(Psi(i,j,:),int);
else
XOUT(i,j,:)=nan*ones([1,p]);
end
end
end
Psitide=XOUT;
Psilow=Psi-Psitide;
T0=17/24;
Fc=1/T0;
dt=1/24;
fN=1/(2*dt); % Nyquist frequency
N=round(2*mean(T0/dt)); % window length
%N=min(N)
B = fir1(N,Fc/fN,'Low'); % filter coefficients
I2 = find(~isnan(Psi));
if ~isempty(I2)
Psilow2= filtfilt(B,1,Psilow);%-nanmean(x));
else
Psilow2= nan*ones(size(Psilow));
end
%interpoler sur grille reguliere si c'est pas fait%%%%%%regular; non rotated! Xr is regular!
%%%%%%%%%
dx=0.05;
dy=dx;
X=LON;
Y=LAT;
xr=min(X(:)):dx:max(X(:));
yr=min(Y(:)):dy:max(Y(:));
[Xr0,Yr0]=meshgrid(xr,yr);
%interpoler sur grille reguliere %%%%%%regular; non rotated! Xr is regular!
% Psir=0*ones([size(Xr0,1), size(Xr0,2), size(Psi,3)]);
for t=1:size(Psilow,3)
Psilow3(:,:,t)=griddata(X,Y,Psilow2(:,:,t),Xr0,Yr0,'linear');
end
[gxpsi,gypsi]=gradient(Psilow3,dx*1000.*111*cosd(mean(Y(:))),dy*1000*111,1);
ugeo=-gypsi;
vgeo= gxpsi;
return
main_RKmA.m 0 → 100755
View file @ f1cc58c8
% 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
% email: tamkpanka.tamtare@uqar.ca
% juin 2018; Last revision:
%------------- BEGIN CODE ------------------------------------------------------------------------------------------------------------------------------------
%*************************************************************************************************************************************************************
% Definition des constantes:
%*************************************************************************************************************************************************************
% h............................Pas de temps pour le transport de particules (1 = 1 tstep
% donc une heure pcq les sorties du modele sont horaires).
% i1...........................Valeur initiale de i pour le positionnement des particu-
% les initiales.
% i2...........................Valeur finale de i pour le positionnement des particules
% initiales.
% k1...........................Valeur initiale de k pour le positionnement des particu-
% les initiales.
% k2...........................Valeur finale de k pour le positionnement des particules
% initiales.
%*************************************************************************************************************************************************************
% Definition des variables:
%*************************************************************************************************************************************************************
% 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
%*************************************************************************************************************************************************************
% 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;
else % Ensemble des cellules
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
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
end
%******************************************************************Fin de la configuration*********************************************************************
%****************************************************************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';
load 'juin12014.mat';
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 = 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=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);
% keyboard
%**********************************************************************************************************************************************************
% 2. Transporter les particules avec Runge-Kutta 4
%**********************************************************************************************************************************************************
tsteps1 = T1:T2-1;
tsteps11=tsteps1+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);
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
else
data = [time,lon,lat];
end
main_RKmB.m 0 → 100755
View file @ f1cc58c8
% 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
% email: tamkpanka.tamtare@uqar.ca
% juin 2018; Last revision:
%------------- BEGIN CODE ------------------------------------------------------------------------------------------------------------------------------------
%*************************************************************************************************************************************************************
% Definition des constantes:
%*************************************************************************************************************************************************************
% h............................Pas de temps pour le transport de particules (1 = 1 tstep
% donc une heure pcq les sorties du modele sont horaires).
% i1...........................Valeur initiale de i pour le positionnement des particu-
% les initiales.
% i2...........................Valeur finale de i pour le positionnement des particules
% initiales.
% k1...........................Valeur initiale de k pour le positionnement des particu-
% les initiales.
% k2...........................Valeur finale de k pour le positionnement des particules
% initiales.
%*************************************************************************************************************************************************************
% Definition des variables:
%*************************************************************************************************************************************************************
% 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
%*************************************************************************************************************************************************************
% 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;
else % Ensemble des cellules
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
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
end
%******************************************************************Fin de la configuration*********************************************************************
%% ************ 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';
load 'juin12014.mat';
dat= datO;
uca=u;
vca=v;
load 'juin22014.mat';
ucb=u;
vcb=v;
%% vent du model GEM (test avec le spot 73)
load 'GEM2014.mat';
uw=Uw;
vw=Vw;
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
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=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);
% keyboard
%****************************************************************************************************************************************************
% Transporter les particules avec Runge-Kutta 4
%****************************************************************************************************************************************************
tsteps1 = T1:T2-1;
tsteps11=tsteps1+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)
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;
[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');
if numel(time) ~= numel(lat)
disp('numel(time) ~= numel(lat)')
whos spot time lat lon
data = [NaN];
return
else
data = [time,lon,lat];
end
main_RKmC.m 0 → 100644
View file @ f1cc58c8
% 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
% email: tamkpanka.tamtare@uqar.ca
% juin 2018; Last revision:
%------------- BEGIN CODE ------------------------------------------------------------------------------------------------------------------------------------
%*************************************************************************************************************************************************************
% Definition des constantes:
%*************************************************************************************************************************************************************
% h............................Pas de temps pour le transport de particules (1 = 1 tstep
% donc une heure pcq les sorties du modele sont horaires).
% i1...........................Valeur initiale de i pour le positionnement des particu-
% les initiales.
% i2...........................Valeur finale de i pour le positionnement des particules
% initiales.
% k1...........................Valeur initiale de k pour le positionnement des particu-
% les initiales.
% k2...........................Valeur finale de k pour le positionnement des particules
% initiales.
%*************************************************************************************************************************************************************
% Definition des variables:
%*************************************************************************************************************************************************************
% k............................Indice de latitude dans la grille du modèle
% iii..........................Indice de longitude dans la grille du modèle
%*************************************************************************************************************************************************************
function [data] = main_RKmC(dat1,dat2,lat1,lon1)
close all
%clear mex
%*************************************************************************************************************************************************************
% 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