Skip to content

GitLab

W
WIM2
Commit d380ec6c authored by Jérémy Baudry's avatar Jérémy Baudry
Browse files
Options

second commit

parent 81dede1c
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 973 additions and 87 deletions
WIM2 0 → 100755
View file @ d380ec6c
File added
nml/parameter.nml
View file @ d380ec6c
......@@ -2,8 +2,8 @@
&waves_parameters
Tm =6
Hs =1
disp =0
Hs =3
disp =1
......@@ -11,9 +11,12 @@ disp =0
&model_parameter
nbin =10
nbin =40
dx =5000
Cfl =0.7
Cfl =1
name_sim ='simulation1'
root = 'output/'
FSD_scheme =1
/
......@@ -31,10 +34,9 @@ gamma_s =3.3
/
&ice_parameters
cice =0.5
cice =0.95
hice =1
D0 =500
gam =1.5
Dmin =20
/
output/Energy_spectrum.dat 0 → 100644
View file @ d380ec6c
This diff is collapsed.
output/advection.m 0 → 100644
View file @ d380ec6c
function E = advection(E,c,dx,dt)
%ADVECTION is a 1d advection code Lax Wendroff scheme with Superbee flux limiter
%
% E = advection(E,c,dx,dt)
%
% INPUTS:
% E = vector of the thing to be advected;
% c = scalar speed;
% dx = the spatial resolution;
% dt = the temporal resolution;
%
% OUTPUT:
% E = vector of the thing to be advected - after advection;
if ~isequal(size(E),size(c))
c = c';
end
f = flux(E,c,dx,dt);
E = E-dt*diffl(f)/dx;
function f = flux(E,c,h,dt)
% Lax-Wendroff with Superbee flux limiting;
theta = diffl(E)./(diffr(E)+3e-14);
phi = limiter(theta);
f = c.*E+c/2.*(1-c*dt/h).*diffr(E).*phi;
end
function y = diffl(x)
y = [x(1);diff(x)];
end
function y = diffr(x)
y = [diff(x);-x(end)];
end
% Superbee
function phi = limiter(r)
phi = max(0,max(min(1,2*r),min(r,2)));
end
end
output/advection.m~ 0 → 100644
View file @ d380ec6c
function E = advection(E,c,dx,dt)
%ADVECTION is a 1d advection code Lax Wendroff scheme with Superbee flux limiter
%
% E = advection(E,c,dx,dt)
%
% INPUTS:
% E = vector of the thing to be advected;
% c = scalar speed;
% dx = the spatial resolution;
% dt = the temporal resolution;
%
% OUTPUT:
% E = vector of the thing to be advected - after advection;
if ~isequal(size(E),size(c))
c = c';
end
f = flux(E,c,dx,dt);
E = E-dt*diffl(f)/dx;
function f = flux(E,c,h,dt)
% Lax-Wendroff with Superbee flux limiting;
theta = diffl(E)./(diffr(E)+3e-14);
phi = limiter(theta);
f = c.*E+c/2.*(1-c*dt/h).*diffr(E).*phi;
end
function y = diffl(x)
y = [0;diff(x)];
end
function y = diffr(x)
y = [diff(x);0];
end
% Superbee
function phi = limiter(r)
phi = max(0,max(min(1,2*r),min(r,2)));
end
end
output/data_treatment.m
View file @ d380ec6c
% clear all;
clear all;
close all;
E=load('output/Energy_spectrum.dat');
subplot(1,2,1)
cmap=jet(size(E,1));
for i=1:size(E,1)
plot(E(i,:),'color',cmap(i,:))
hold on
x=ncread('simulation1.nc','x_axis');
t=ncread('simulation1.nc','time');
om=ncread('simulation1.nc','omega');
spectre=ncread('simulation1.nc','Spectrum');
Dave=ncread('simulation1.nc','Dave');
Dmax=ncread('simulation1.nc','Dmax');
f=om/(2*pi);
E=reshape(spectre(end,end,:),length(om),1);
Ei=reshape(spectre(30,:,:),length(x),length(om));
E1=reshape(spectre(1,1,:),length(om),1);
figure(1)
cmap=rand(length(t),3);
w = waitforbuttonpress;
for i=1:length(t)
figure(1)
sp=reshape(spectre(i,:,:),length(x),length(om));
h=mesh(f,x,sp);
axis([min(f) max(f) min(x) max(x) 0 max(E1)])
xlabel('Frequency [s^{-1}]')
zlabel('Energy')
ylabel('x [km]')
pause(0.1)
end
FSD=load('output/floe_size.dat');
subplot(1,2,2)
plot(FSD(:,1))
for i=1:length(om)
EE(i)=sum(Ei(:,i));
end
figure
plot(E1)
hold on
plot(FSD(:,2),'r')
plot(EE,'r')
for i=1:length(t)
E2=reshape(spectre(i,:,:),length(x),length(om));
m0(i)=sum(sum(E2))/sum(E1);
end
figure
plot(t,m0)
figure
plot(x,Dmax,'color','r')
hold on
plot(x,Dave,'--b')
grid on
xlabel('x [km]')
ylabel('Floe size [m]')
legend('Dmax','<D>')
\ No newline at end of file
output/data_treatment.m~ 0 → 100644
View file @ d380ec6c
clear all;
close all;
x=ncread('simulation1.nc','x_axis');
t=ncread('simulation1.nc','time');
om=ncread('simulation1.nc','omega');
spectre=ncread('simulation1.nc','Spectrum');
Dave=ncread('simulation1.nc','Dave');
Dmax=ncread('simulation1.nc','Dmax');
f=om/(2*pi);
E=reshape(spectre(end,end,:),length(om),1);
Ei=reshape(spectre(30,:,:),length(x),length(om));
E1=reshape(spectre(1,1,:),length(om),1);
figure(1)
cmap=rand(length(t),3);
w = waitforbuttonpress;
for i=1:length(t)
figure(1)
sp=reshape(spectre(i,:,:),length(x),length(om));
h=mesh(f,x,sp);
axis([min(f) max(f) min(x) max(x) 0 max(E1)])
xlabel('Frequency [s^{-1}]')
zlabel('Energy')
ylabel('x [km]')
pause(0.1)
end
for i=1:length(om)
EE(i)=sum(Ei(:,i));
end
figure
plot(E1)
hold on
plot(EE,'r')
for
% figure
%
% subplot(1,2,1)
%
% plot(T,E,'color','r')
% hold on
% plot(T,Ei,'--b')
% grid on
% xlabel('Frequency [s^{-1}]')
% ylabel('Energy')
%
%
% subplot(1,2,2)
%
% plot(x,Dmax,'color','r')
% hold on
% plot(x,Dave,'--b')
% grid on
%
% xlabel('x [km]')
% ylabel('Floe size [m]')
%
% legend('Dmax','<D>')
\ No newline at end of file
output/ferret.jnl 0 → 100644
View file @ d380ec6c
! NOAA/PMEL TMAP
! FERRET v6.93
! Linux 2.6.32-504.el6.x86_64 64-bit - 11/13/14
! 1-Oct-15 15:51
use simulation1
use simulation1.nc
show data
use simulation1.nc
exit
output/ferret.jnl.~1~ 0 → 100644
View file @ d380ec6c
! NOAA/PMEL TMAP
! FERRET v6.93
! Linux 2.6.32-504.el6.x86_64 64-bit - 11/13/14
! 30-Sep-15 15:02
use simulation1.nc
use simulation1
quit
output/ferret.jnl.~2~ 0 → 100644
View file @ d380ec6c
! NOAA/PMEL TMAP
! FERRET v6.93
! Linux 2.6.32-504.el6.x86_64 64-bit - 11/13/14
! 30-Sep-15 15:05
use simulation1.nc
show data
show axis/all
show data
plot/I=1/J=1 SPECTRUM
output/floe_size.dat 0 → 100644
View file @ d380ec6c
0.0000000000000000 0.0000000000000000
500.00000000000000 500.00000000000000
500.00000000000000 500.00000000000000
500.00000000000000 500.00000000000000
500.00000000000000 500.00000000000000
500.00000000000000 500.00000000000000
500.00000000000000 500.00000000000000
500.00000000000000 500.00000000000000
500.00000000000000 500.00000000000000
500.00000000000000 500.00000000000000
output/jonswap.m 0 → 100755
View file @ d380ec6c
function [E] = JONSWAP(Tm,Hs,om)
% Formule http://hmf.enseeiht.fr/travaux/CD0910/bei/beiere/groupe4/node/59
gam = 3.3;
alp = 0.0624/(0.23+0.0336*gam-0.185/(1.9+gam));
fp = 1/Tm;
f = om/(2*pi);
sig = zeros(length(om),1);
sig(f<=fp) = 0.07;
sig(f>fp) = 0.09;
E = alp*Hs^2*fp^4*f.^-5.*exp(-5/4*(fp./f).^4).*gam.^(exp(-(f-fp).^2./(2.*sig.^2.*fp^2)));
end
output/test_adv 0 → 100644
View file @ d380ec6c
clear all;
close all;
% Parameters
g = 9.81; % gravitational acceleration
dx = 5000; % spatial resolution
Hs=1;
Tm=6;
% Waves
fmin = 1/20; % minimum wave frequency
fmax = 1/2.5; % maximum wave frequency
om1 = 2*pi*fmin; % minimum wave radial frequency
om2 = 2*pi*fmax; % maximum wave radial fequency
nw = 61; % number of frequency bins
dw = (om2-om1)/(nw-1); % integral interval for wave radial frequencies
om = om1+(0:nw-1)'*dw; % wave radial frequencies vector
T = 2*pi./om; % wave periods vector
wlng = g.*T.^2./(2.*pi); % wavelength as a function only of wave period
cp = sqrt(g.*wlng./(2.*pi)); % phase speed
cg = cp./2; % group speed
cgmax = max(cg); % group speed maximum
cg(:) = cgmax; % no dispersion = all group speed are the same (maximum)
dt = dx/cgmax; % time interval (temporal resolution)
Ei = jonswap(Tm,Hs,om); % JONSWAP spectrum
\ No newline at end of file
output/test_adv.m 0 → 100644
View file @ d380ec6c
clear all;
close all;
% Parameters
g = 9.81; % gravitational acceleration
dx = 500; % spatial resolution
nbin=80;
x=linspace(0,dx*nbin/1000,nbin);
Hs=1;
Tm=6;
% Waves
fmin = 1/20; % minimum wave frequency
fmax = 1/2.5; % maximum wave frequency
om1 = 2*pi*fmin; % minimum wave radial frequency
om2 = 2*pi*fmax; % maximum wave radial fequency
nw = 61; % number of frequency bins
dw = (om2-om1)/(nw-1); % integral interval for wave radial frequencies
om = om1+(0:nw-1)'*dw; % wave radial frequencies vector
T = 2*pi./om; % wave periods vector
wlng = g.*T.^2./(2.*pi); % wavelength as a function only of wave period
cp = sqrt(g.*wlng./(2.*pi)); % phase speed
cg = cp./2; % group speed
cgmax = max(cg); % group speed maximum
cg(:) = cgmax; % no dispersion = all group speed are the same (maximum)
CN=0.7;
dt = CN*dx/cgmax; % time interval (temporal resolution)
nsteps=ceil(nbin/CN);
Ei = jonswap(Tm,Hs,om); % JONSWAP spectrum
E=zeros(nsteps,nbin,length(om));
E(1,1,:)=Ei;
EE=reshape(E(1,:,:),nbin,length(om));
for n=2:nsteps
for w=1:nw % Advection loop over each frequency
EE(:,w) = advection(EE(:,w),cg(w,:),dx,dt); % SEE advection routine
end
E(n,:,:)=reshape(EE,1,nbin,length(om));
end
figure(1)
rr = waitforbuttonpress;
for n=1:nsteps
figure(1)
E3=reshape(E(n,:,:),nbin,length(om));
h=mesh(om,x,E3);
axis([om1 om2 min(x) max(x) 0 max(Ei)])
xlabel('Frequency [s^{-1}]')
zlabel('Energy')
ylabel('x [km]')
pause(0.1)
end
for i=1:nsteps
E2=reshape(E(i,:,:),nbin,length(om));
m0(i)=sum(sum(E2))/sum(Ei);
end
figure
plot(m0)
output/test_adv.m~ 0 → 100644
View file @ d380ec6c
clear all;
close all;
% Parameters
g = 9.81; % gravitational acceleration
dx = 500; % spatial resolution
nbin=80;
x=linspace(0,dx*nbin/1000,nbin);
Hs=1;
Tm=6;
% Waves
fmin = 1/20; % minimum wave frequency
fmax = 1/2.5; % maximum wave frequency
om1 = 2*pi*fmin; % minimum wave radial frequency
om2 = 2*pi*fmax; % maximum wave radial fequency
nw = 61; % number of frequency bins
dw = (om2-om1)/(nw-1); % integral interval for wave radial frequencies
om = om1+(0:nw-1)'*dw; % wave radial frequencies vector
T = 2*pi./om; % wave periods vector
wlng = g.*T.^2./(2.*pi); % wavelength as a function only of wave period
cp = sqrt(g.*wlng./(2.*pi)); % phase speed
cg = cp./2; % group speed
cgmax = max(cg); % group speed maximum
cg(:) = cgmax; % no dispersion = all group speed are the same (maximum)
CN=0.5;
dt = CN*dx/cgmax; % time interval (temporal resolution)
nsteps=ceil(nbin/CN);
Ei = jonswap(Tm,Hs,om); % JONSWAP spectrum
E=zeros(nsteps,nbin,length(om));
E(1,1,:)=Ei;
EE=reshape(E(1,:,:),nbin,length(om));
figure(1)
rr = waitforbuttonpress;
for n=2:nsteps
for w=1:nw % Advection loop over each frequency
EE(:,w) = advection(EE(:,w),cg(w,:),dx,dt); % SEE advection routine
end
% Incident wave spectrum
%EE(1,:) = Ei;
figure(1)
h=mesh(om,x,EE);
axis([om1 om2 min(x) max(x) 0 max(Ei)])
xlabel('Frequency [s^{-1}]')
zlabel('Energy')
ylabel('x [km]')
pause(0.1)
E(n,:,:)=reshape(EE,1,nbin,length(om));
end
for i=1:nsteps
E2=reshape(E(i,:,:),nbin,length(om));
m0(i)=sum(sum(E2))/sum(Ei);
end
figure
plot(m0)
output/wim.m 0 → 100755
View file @ d380ec6c
function [T,cice,Hsig,Ei,E,S_win,S_wcp,S_ice,Dmax,Dave,om_do] = wim(C,h,D,U10,Tm,Hs)
% WIM - Models the propagation of a wave spectrum along a 1D
% ice-covered transect.
%
% Syntax: [T,cice,Hsig,Ei,E,S_win,S_wcp,S_ice] = wim_final(C,h,D,U10,Tm,Hs)
%
% Inputs:
% C : ice concentration. For a homogeneous concentration, C is scalar
% in tenths. Otherwise, C is a vector of the size of the spatial
% grid, in fraction of 1 (0 < C < 1). In this version of WIM, if
% C is homogeneous, the transect is 5km long with spatial
% resolution of 500m.
% h : ice thickness in meters. h is a scalar.
% D : average floe diameter in meters. D is a scalar.
% U10 : wind speed at 10m elevation in m/sec. U10 is a scalar.
% Tm : mean wave period to build the spectrum, in seconds. Tm is a scalar.
% Hs : significant wave height, in meters. Hs is a scalar.
%
% Outputs:
% T : wave period in seconds. T is a vector the size of the frequency range.
% cice : ice concentration vector in fraction of 1. cice has the size of
% the spatial grid.If C is a vector, then cice = C.
% Hsig : final wave heigh in each cell in meters. Hsig has the size of
% the spatial grid.
% Ei : initial spectrum in m2/hz.
% E : final wave spectrum in each cell [m2/Hz].
% S_win: final wind source term in each cell in m2/Hz/sec.
% S_wcp: final white-capping source term in each cell in m2/Hz/sec.
% S_ice: final ice source term in each cell in m2/Hz/sec.
% Parameters
g = 9.81; % gravitational acceleration
dx = 5000; % spatial resolution
% Ice conditions
if length(C) == 1 % Homogeneous ice concentration case
nx = 10; % 10 cells of 500m -> 5km long transect
cice = zeros(1,nx)+C; % homogeneous ice concentration vector
hice = zeros(1,nx)+h; % homogeneous ice thickness vector
Dmax = zeros(1,nx)+D; % homogeneous floe diameter vector
Dave = zeros(1,nx)+D;
else % Variable ice concentration case
nx = length(C); % number of cells
cice = C; % ice concentration vector
hice = zeros(1,nx); % ice thickness vector size allocation
hice(find(C)) = h; % ice thickness vector
Dmax = zeros(1,nx); % floes diameter vector size allocation
Dmax(find(C)) = D; % floes diameter vector
end
% Waves
fmin = 1/20; % minimum wave frequency
fmax = 1/2.5; % maximum wave frequency
om1 = 2*pi*fmin; % minimum wave radial frequency
om2 = 2*pi*fmax; % maximum wave radial fequency
nw = 61; % number of frequency bins
dw = (om2-om1)/(nw-1); % integral interval for wave radial frequencies
om = om1+(0:nw-1)'*dw; % wave radial frequencies vector
T = 2*pi./om; % wave periods vector
Ei = jonswap(Tm,Hs,om); % JONSWAP spectrum
wlng = g.*T.^2./(2.*pi); % wavelength as a function only of wave period
cp = sqrt(g.*wlng./(2.*pi)); % phase speed
cg = cp./2; % group speed
cgmax = max(cg); % group speed maximum
cg(:) = cgmax; % no dispersion = all group speed are the same (maximum)
% Wind
U10 = repmat(U10,nx,1); % wind speed at 10m elevation vector
% Temporal grid
dt = dx/cgmax; % time interval (temporal resolution)
nsteps = nx+1; % number of time steps (=nx+1 because advection is
% done from one cell at each time step)
time = 0:dt:nsteps*dt; % time vector
nt = length(time); % time vector size
% Memory preallocation
E = zeros(nx,nw); % wave spectrum
Swin = zeros(nx,nw); % wind source term
Swcp = zeros(nx,nw); % white-capping source term
Sice = zeros(nx,nw); % ice dissipation source term
S_win = zeros(nx,nw); % wind source term weighted by open water fraction (for outputs)
S_ice = zeros(nx,nw); % ice dissipation source terme weighted by ice fraction
S_wcp = zeros(nx,nw); % white-capping source term weighted by open water fraction
Hsig = zeros(1,nx); % significant wave height
% Action density limiter
Slim = action_density_limiter(om,cp); % SEE action_density_limiter routine
for n=1:nt % Time loop
% Advection.
% The advection is done by solving Dt(S) = 0 using the
% Lax-Wendroff scheme with Superbee flux limiting and a Neumann
% boundary condition. The advection is performed over the whole domain
% in one step on an unattenuated intermediate spectrum.
for w=1:nw % Advection loop over each frequency
E(:,w) = advection(E(:,w),cg(w,:),dx,dt); % SEE advection routine
end
% Incident wave spectrum
E(1,:) = Ei; % The initial wave spectrum is forced in the first cell
% Processes for each spatial cell with ice (generation by wind, dissipation
% by white-capping and attenuation by ice)
for i=2:nx % Spatial loop
% % Generation by wind
% Swin(i,:) = wind_gen(U10(i),E(i,:),om,cp); % SEE wind_gen routine
% Swin(i,:) = min(Swin(i,:),Slim'); % action density limiter
% E(i,:) = E(i,:) + Swin(i,:)*dt*(1-cice(i)); % the wave spectrum is updated (explicit scheme)
% S_win(i,:) = Swin(i,:)*(1-cice(i)); % effective wind source term for outputs
% % Dissipation by white-capping
% Swcp(i,:) = white_cap(E(i,:),om,cp); % SEE white-cap routine
% Swcp(i,:) = max(Swcp(i,:),-Slim'); % action density limiter
% E(i,:) = E(i,:) + Swcp(i,:)*dt*(1-cice(i)); % the wave spectrum is updated (semi-implicit scheme)
% S_wcp(i,:) = Swcp(i,:)*(1-cice(i)); % effective white-capping source term for outputs
% Attenuation by ice
if cice(i)>0
Sice(i,:) = ice_att(E(i,:),T,cg,cice(i),hice(i),Dave(i),dt*cice(i)); % SEE ice_att routine
E(i,:) = E(i,:) + Sice(i,:)*dt; % the wave spectrum is updated (implicit scheme)
else
Sice(i,:) = 0; % if there is no ice, there is no attenuation Sherlock !
end
S_ice(i,:) = Sice(i,:)*cice(i); % effective ice dissipation source term for outputs
[Dmax(i),om_do(i)]=floe_breaking(E(i,:),h,om,Dmax(i),dw);
Dave(i)=FSD(Dmax(i),D);
% Wave spectrum statistics
E(E<0) = 0; % avoiding negative energy
m0 = trapz(om,E(i,:)); % total energy
Hsig(i) = 4*sqrt(m0); % significant wave height
end
end
end
src/advection.f90
View file @ d380ec6c
......@@ -4,35 +4,59 @@ use parameters
implicit none
double precision, allocatable :: beta_1(:),beta_2(:),beta_3(:)
allocate(beta_1(nfreq))
allocate(beta_2(nfreq))
allocate(beta_3(nfreq))
beta_1=CN/2*(CN+1)
beta_2=1-CN**2