Commit 5e256d18 by Jérémy Baudry

### nettoyage!

parent d380ec6c
WIM2 deleted 100755 → 0
File deleted
This diff is collapsed.
File deleted
 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
 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
 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','') \ No newline at end of file
 ! 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
 ! 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
 ! 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
 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
 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
File deleted
 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
 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)
 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)
File deleted
File deleted
File deleted
File deleted
File deleted
File deleted
 COMPILER= ifort OPTION= -O3 NETCDFinc= -I/usr/local/include NETCDFLIB= -L/usr/local/lib -lnetcdf OBJ= *.o SRC= *.f90 MOD= parameters.mod EXEC= WIM2 WIM2: \$(OBJ) \$(COMPILER) \$(OPTION) \$(NETCDFinc) -o \$(EXEC) \$(OBJ) \$(NETCDFLIB) mv WIM2 ../ \$(MOD): parameters.f90 \$(COMPILER) -c parameters.f90 \$(OBJ): \$(SRC) \$(MOD) \$(COMPILER) \$(NETCDFinc) -c \$(SRC) \$(NETCDFLIB) clean: rm -f *.o *.mod mrproper: rm -f *.o *.mod ../WIM2
This diff is collapsed.
File deleted
File deleted
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!