Skip to content

GitLab

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

new version

parent 5e256d18
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 528 additions and 120 deletions
nml/parameter.nml
View file @ 4508a952
......@@ -2,8 +2,8 @@
&waves_parameters
Tm =6
Hs =3
disp =1
Hs =1.5
disp =0
......@@ -11,8 +11,8 @@ disp =1
&model_parameter
nbin =40
dx =5000
nbin =100
dx =1000
Cfl =1
name_sim ='simulation1'
root = 'output/'
......@@ -22,8 +22,8 @@ FSD_scheme =1
&spectrum_parameters
nfreq =60
init_spec =1
nfreq =200
Tmin =2.5
Tmax =20
alpha_s =0.2044
......@@ -34,8 +34,12 @@ gamma_s =3.3
/
&ice_parameters
cice =0.95
hice =1
X_ice =5000
cice =0.8
ice_thick =0
hice =2
hmax =4
Xh =100000
D0 =500
gam =1.5
Dmin =20
......
output/data_treatment.m
View file @ 4508a952
......@@ -8,67 +8,64 @@ 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 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
hice=ncread('simulation1.nc','Ice thickness');
c=ncread('simulation1.nc','Ice concentration');
FSD=ncread('simulation1.nc','Floe size distribution');
floesize=ncread('simulation1.nc','floe size');
% 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
%
%
%
%
% figure
%
% plot(x,hice,'color','r')
% hold on
% plot(x,c,'--b')
% grid on
%
% xlabel('x [km]')
%
% legend('ice thickness','ice concentration')
%
%
%
% 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/floe_size.dat 0 → 100644
View file @ 4508a952
This diff is collapsed.
src/attenuation.f90
View file @ 4508a952
......@@ -36,15 +36,13 @@ p2 = q21*T**4 + q22*T**3 + q23*T**2 + q24*T + q25
p3 = q31*T**4 + q32*T**3 + q33*T**2 + q34*T + q35
alpha(i,1:nfreq)=-1*(p1*h(i)**2 + p2*h(i) + p3)
do j=1,nfreq
if (alpha(i,j).lt.0)then
alpha(i,j)=0
alpha(i,j)=0d0
end if
end do
att=C_ice(i)*alpha(i,1:nfreq)/Dave(i)
att=C_ice(i)*alpha(i,1:nfreq)/(Dave(i)+3e-14)
S_ice(i,1:nfreq)=-att*E(n,i,1:nfreq)
......
src/break_horvat.f90 0 → 100644
View file @ 4508a952
subroutine break_horvat
use parameters
implicit none
integer::fl,iii,jjj,bin_new_floe
double precision :: rho=1025
double precision :: strain_crit=3e-5,new_floe
double precision, allocatable :: Amp(:),P_wa(:)
double precision, allocatable :: strain(:),P_f(:,:),coef(:)
allocate(Amp(nfreq))
allocate(P_wa(nfreq))
allocate(strain(nfreq))
allocate(P_f(nbcat,nfreq))
allocate(coef(nbcat))
P_f=0
!transform energy spectrum into amplitude spectrum
Amp(1:nfreq)=sqrt(2*E(n,i,1:nfreq)*domega/(rho*g))
!calcul the strain at each wavelength
strain(1:nfreq)=Amp(1:nfreq)*h*2*pi**2/wl(1:nfreq)**2
!find the probability of occurence of these amplitude using rayleigh
!distribution
P_wa(1:nfreq)=2*Amp(1:nfreq)*exp(-Amp(1:nfreq)**2/Hs**2)/Hs**2
!Compute breaking probability
fl=0
do iii=1,nedge-1
fl=fl+1
do jjj=1,nfreq
if (strain(jjj).gt.strain_crit .and. wl(jjj).lt.floe_cat(iii)) then
P_f(fl,jjj)=P_wa(jjj)
else
P_f(fl,jjj)=0
end if
end do
!coef(fl)=sum(P_f(fl,1:nfreq))
!P_f(fl,1:nfreq)=P_f(fl,1:nfreq)/coef(fl)
end do
! do the fragmentation
fl=0
do iii=1,nedge-1
fl=fl+1
do jjj=1,nfreq
if (strain(jjj).gt.strain_crit .and. wl(jjj).lt.floe_cat(iii)) then
new_floe=0.5*wl(jjj)
bin_new_floe=nedge-ceiling((new_floe-minfloe)/res)
FSD(n,i,bin_new_floe)=FSD(n,i,bin_new_floe)+P_f(fl,jjj)*FSD(n,i,fl)
FSD(n,i,fl)=FSD(n,i,fl)-P_f(fl,jjj)*FSD(n,i,fl)
end if
end do
end do
if(i.eq.2 .and. n.eq.2) then
do fl=1,nbcat
write(11,*)P_f(fl,1:nfreq)
end do
end if
end subroutine break_horvat
src/floe_breaking.f90
View file @ 4508a952
......@@ -5,7 +5,7 @@ use parameters
implicit none
double precision, allocatable :: WN_ice(:),fnt(:),k(:)
double precision ::rho=1025,d,F,Y=5.5,v=0.3
integer :: acc=1000
integer :: acc=10000
double precision, allocatable :: W(:)
double precision, allocatable :: fnt_d(:)
......@@ -28,12 +28,12 @@ implicit none
allocate(S(nfreq))
F=Y*h(i)**3/12*(1-v**2)
F=Y*h(i)**3/12d0*(1-v**2)
d=0.9*h(i)
jj=1
do ii=0,acc
k(jj)=real(ii)*2/real(acc)
k(jj)=real(ii)*10d0/real(acc)
jj=jj+1
end do
......@@ -46,12 +46,12 @@ implicit none
W=g*WN_ice/omega**2
S=(h(i)/2)*WN_ice**2*W
S=(0.5*h(i))*WN_ice**2*W
m0_s=sum(E(n,i,1:nfreq)*S**2)*domega
SS=2*sqrt(m0_s)
if (SS.gt.RS) then
m0_a=sum(E(n,i,1:nfreq)*W**2)*domega
m2_a=sum(omega**2*E(n,i,1:nfreq)*W**2)*domega
Tw=2*pi*sqrt(m0_a/m2_a)
......@@ -62,8 +62,8 @@ implicit none
wl_w=2*pi/WN_ice_d
Dmax(i)=max(Dmin,min(wl_w/2,Dmax(i)))
Dmax(i)=max(Dmin,min(0.5*wl_w,Dmax(i)))
end if
......
src/fsd_build.f90
View file @ 4508a952
......@@ -3,16 +3,43 @@ use parameters
implicit none
double precision :: coeff
double precision, allocatable :: ND(:),NN(:)
integer :: M,mm
if (Dmax(i).eq.Dmin) then
Dave(i)=Dmin
elseif (Dmax(i).eq.D0) then
Dave(i)=Dmax(i)
if(C_ice(i).eq.0)then
Dave(i)=0
elseif (Dmax(i).eq.D0) then
Dave(i)=Dmax(i)
elseif (Dmax(i).eq.Dmin) then
Dave(i)=Dmin
elseif (FSD_scheme.eq.1) then
coeff=1/((1/(1-gam))*(Dmax(i)**(1-gam)-Dmin**(1-gam)))
Dave(i)=coeff*(1/(2-gam))*(Dmax(i)**(2-gam)-Dmin**(2-gam))
else
M=floor(log(Dmax(i)/Dmin)/log(psi))
allocate(ND(M+1))
allocate(NN(M+1))
jj=1
do mm=0,M-1
NN(jj)=((1-ff)*(psi**2*ff)**mm)*(Dmax(i)/psi**mm)
ND(jj)=(1-ff)*(psi**2*ff)**mm
jj=jj+1
end do
Dave(i)=(sum(NN)+((psi**2*ff)**M)*(Dmax(i)/psi**M))/(sum(ND)+(psi**2*ff)**M)
deallocate(ND)
deallocate(NN)
end if
end if
end subroutine fsd_build
src/graph_dislin.f90 0 → 100644
View file @ 4508a952
subroutine graph_dislin
use parameters
use dislin
IMPLICIT NONE
double precision, allocatable :: axe(:),Ep(:,:)
real ::par1,par2,par3,par4,par5,par6,par7,par8,par9,par10,par11,par12
allocate(axe(nbin))
allocate(Ep(nbin,nfreq))
axe(1)=0
do i=1,nbin-1
axe(i+1)=i*dx/1000
end do
!_____________________________________maximum floe size________________
!CALL SETPAG('DA4P')
!CALL METAFL('CONS')
!CALL SCRMOD('REVERS')
!CALL DISINI()
!CALL PAGERA()
!CALL SETRGB (0.5,0.2,0.4)
!CALL GRAF (real(axe(1),4),real(axe(nbin),4),real(0,4),real(10*dx/1000,4),real(0,4),real(maxval(Dmax),4),real(0,4),real(50,4))
!call NAME('x axis [km]','X')
!call NAME('Dmax [m]','Y')
!CALL TITLIN ('Maximum Floe size in the MIZ',1)
!call TITLE
! call COLOR('GRAY')
! call GRID(10,5)
! call COLOR('RED')
! CALL CURVE (real(axe,4),real(Dmax,4),nbin)
! call ENDGRF
!call disfin
!________________________SPECTRUM_______________________________________
do i=1,nbin
do ii=1,nfreq
Ep(i,ii)=real(E(15,i,ii),4)
end do
end do
par5=real(omega(1),4)
par6=real(omega(nfreq),4)
par7=real(omega(1),4)
par8=real(1,4)
par1=real(axe(1),4)
par2=real(axe(nbin),4)
par3=real(0,4)
par4=real(10*dx/1000,4)
par9=real(0.,4)
par10=real(maxval(Ei),4)
par11=real(0.,4)
par12=real(maxval(Ei)/10,4)
!write(*,*)par1,par2,par3,par4,par5,par6,par7,par8,par9,par10,par11,par12
CALL SETPAG('DA4P')
CALL METAFL('CONS')
!CALL SCRMOD('REVERS')
CALL DISINI()
CALL PAGERA()
CALL VIEW3D(-5.,-5.,4.,'ABS')
do n=1,nsteps
do i=1,nbin
do ii=1,nfreq
Ep(i,ii)=real(E(n,i,ii),4)
end do
end do
call ERASE
call COLOR('WHITE')
CALL GRAF3D (par1,par2,par3,par4,par5,par6,par7,par8,par9,par10,par11,par12)
CALL NAME('X-axis [km]','X')
CALL NAME('Frequency','Y')
CALL NAME('Energy','Z')
call COLOR('BLUE')
CALL SURMAT (real(Ep,4),nbin,nfreq,1,1)
call ENDGRF
call system('sleep 0.1')
end do
call disfin
end subroutine graph_dislin
src/initialization.f90
View file @ 4508a952
......@@ -4,6 +4,7 @@ use parameters
!local parameters
implicit none
double precision, allocatable ::Gf(:),PM(:)
integer ::X1
allocate(Gf(nfreq))
allocate(PM(nfreq))
......@@ -12,7 +13,8 @@ implicit none
!_________________________INITIAL SPECTRUM_____________________________
E(1:nsteps,1:nbin,1:nfreq)=0d0
if(init_spec.eq.1) then
!build JONSWAP spectrum
do i=1,nfreq
......@@ -27,6 +29,11 @@ implicit none
PM=alpha_s*Hs**2*(freq_s**4/freq**5)*exp(-beta_s*(freq_s/freq)**4)
Ei=Gf*PM
else
Ei=(1.25*Hs**2*(1/freq)**5)/(8*pi*Tm**4)*exp(-1.25*((1/freq)/Tm)**4)
end if
......@@ -34,4 +41,35 @@ implicit none
!_______________________________________________________________________
!_______________________ICE_TRANSECT____________________________________
X1=floor(X_ice/dx)
C_ice(1:X1)=0
C_ice(X1:nbin)=cice
Dave(1:X1)=0
Dmax(1:X1)=0
Dave(X1:nbin)=D0
Dmax(X1:nbin)=D0
if (ice_thick.eq.0) then
h(X1:nbin)=hice
h(1:X1)=0d0
else
h(1:X1)=0d0
do jj=X1,nbin
h(jj)=hmax*(0.1+0.9*(1-exp(-((real(jj)*dx-X_ice)/Xh))))
end do
end if
!_____________________FSD_____________________________________________
FSD(1,1:nbin,1)=1d0
end subroutine initialization
src/main.f90
View file @ 4508a952
......@@ -7,7 +7,7 @@ PROGRAM WIM2
call initialization
!______________________________________________________
......@@ -15,6 +15,12 @@ PROGRAM WIM2
spectrum=trim(root)//'Energy_spectrum.dat'
floe_size=trim(root)//'floe_size.dat'
namefile=trim(root)//trim(name_sim)//'.nc'
open(10,file=spectrum)
open(11,file=floe_size)
!_______________________TIME LOOP______________________
do n=2,nsteps
......@@ -23,23 +29,18 @@ PROGRAM WIM2
call advection
end do
do i=2,nbin
do i=1,nbin
call attenuation
call floe_breaking
call fsd_build
call break_horvat
!call floe_breaking
!call fsd_build
end do
end do
!______________________OUTPUTS_________________________
spectrum=trim(root)//'Energy_spectrum.dat'
floe_size=trim(root)//'floe_size.dat'
namefile=trim(root)//trim(name_sim)//'.nc'
open(10,file=spectrum)
open(11,file=floe_size)
call write_output
!call graph_dislin
close(10)
close(11)
......
src/makefile
View file @ 4508a952
......@@ -5,7 +5,7 @@ OPTION= -O3
NETCDFinc=-I/usr/include
NETCDFLIB=-L/usr/lib64 -lnetcdff
NETCDFMOD=-I/usr/lib64/gfortran/modules
DISLIN=-I/usr/local/dislin/gf -ldislin
OBJ= *.o
SRC= *.f90
......@@ -17,7 +17,7 @@ EXEC= WIM2
WIM2: $(OBJ)
$(COMPILER) $(OPTION) $(NETCDFinc) $(NETCDFMOD) -o $(EXEC) $(OBJ) $(NETCDFLIB)
$(COMPILER) $(OPTION) $(NETCDFinc) $(NETCDFMOD) -o $(EXEC) $(OBJ) $(NETCDFLIB) $(DISLIN)
mv WIM2 ../
......@@ -26,7 +26,7 @@ $(MOD): parameters.f90
$(OBJ): $(SRC) $(MOD)
$(COMPILER) $(NETCDFinc) $(NETCDFMOD) -c $(SRC) $(NETCDFLIB)
$(COMPILER) $(NETCDFinc) $(NETCDFMOD) -c $(SRC) $(NETCDFLIB) $(DISLIN)
clean:
......
src/parameters.f90
View file @ 4508a952
......@@ -18,9 +18,11 @@ implicit none
double precision :: Tmin=2.5,Tmax=20,Tm=6,Hs=1
double precision :: gamma_s=3.3,freq_s,domega
double precision, allocatable :: x_axis(:),time(:)
integer :: init_spec=1
!____________waves_________________________________
double precision, allocatable ::wl(:),Cp(:),Cg(:),CN(:)
integer ::disp=0
double precision ::dwl
!____________grid parameters________________________
double precision :: dt,dx=500
integer :: nsteps
......@@ -37,6 +39,21 @@ implicit none
double precision :: D0=200
double precision :: gam=2
double precision :: Dmin=20
double precision :: psi=2
double precision :: ff=0.5
double precision :: Xh=60000
double precision :: X_ice=50000
integer :: ice_thick=1
double precision :: hmax=4
!_________________FSD______________________________
double precision :: minfloe
double precision :: maxfloe
double precision :: res
integer :: nbcat
integer :: nedge
double precision, allocatable :: FSD(:,:,:)