Commit 68586e033d5957e0df6a5f7c53b51c62d07e88c2

Authored by Jérémy Baudry
1 parent a4dc26f5
Exists in master

new release

WIM 0 → 100755
No preview for this file type
batch/.RUN_BATCH_Ec.sh.swp 0 → 100644
No preview for this file type
batch/RUN_BATCH_Ec.sh 0 → 100755
... ... @@ -0,0 +1,35 @@
  1 +#!/bin/bash
  2 +
  3 +
  4 +
  5 +strain="3e-5 4e-5 5e-5 6e-5 7e-5 8e-5 9e-5 1e-4 2e-4 3e-4"
  6 +pc="0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55"
  7 +
  8 +let "count=0"
  9 +for a in $strain
  10 +do
  11 +for b in $pc
  12 +do
  13 + name_exp=pc_ec$count
  14 + cat parameter.txt | sed \
  15 + -e "s/nmebatch/$name_exp/g" \
  16 + -e "s/probcritbatch/$b/g" \
  17 + -e "s/strainbatch/$a/g" \
  18 + -e "s/Hsbatch/1/g" \
  19 + -e "s/hbatch/3/g" \
  20 + -e "s/Tmbatch/8/g" \
  21 + > parameter.nml
  22 +
  23 + mv parameter.nml ../nml/
  24 +
  25 +
  26 + qsub ./lanceur.sh
  27 +
  28 + sleep 30s
  29 +
  30 +
  31 + let "count=count+1"
  32 +done
  33 +done
  34 +
  35 +
... ...
batch/RUN_BATCH_Hs.sh 0 → 100755
... ... @@ -0,0 +1,33 @@
  1 +#!/bin/bash
  2 +
  3 +
  4 +name_exp=('hs1' 'hs2' 'hs3' 'hs4' 'hs5'
  5 + 'hs6' 'hs7' 'hs8' 'hs9' 'hs10')
  6 +
  7 +hs="0.5 1 1.5 2 2.5 3 3.5 4 4.5 5"
  8 +
  9 +let "count=0"
  10 +for a in $hs
  11 +do
  12 + cat parameter.txt | sed \
  13 + -e "s/nmebatch/${name_exp[$count]}/g" \
  14 + -e "s/probcritbatch/0.15/g" \
  15 + -e "s/strainbatch/3e-5/g" \
  16 + -e "s/Hsbatch/$a/g" \
  17 + -e "s/hbatch/4/g" \
  18 + -e "s/Tmbatch/8/g" \
  19 + > parameter.nml
  20 +
  21 + mv parameter.nml ../nml/
  22 +
  23 +
  24 + qsub ./lanceur.sh
  25 +
  26 + sleep 5s
  27 +
  28 +
  29 + let "count=count+1"
  30 +done
  31 +
  32 +
  33 +
... ...
batch/RUN_BATCH_Pc.sh 0 → 100755
... ... @@ -0,0 +1,29 @@
  1 +#!/bin/bash
  2 +
  3 +
  4 +name_exp=('test1' 'test2' 'test3' 'test4')
  5 +
  6 +prob="0.15 0.20 0.30 0.50"
  7 +
  8 +let "count=0"
  9 +for a in $prob
  10 +do
  11 + cat parameter.txt | sed \
  12 + -e "s/nmebatch/${name_exp[$count]}/g" \
  13 + -e "s/probcritbatch/$a/g" \
  14 + -e "s/strainbatch/3e-5/g" \
  15 + > parameter.nml
  16 +
  17 + mv parameter.nml ../nml/
  18 +
  19 +
  20 + qsub ./lanceur.sh
  21 +
  22 + sleep 5s
  23 +
  24 +
  25 + let "count=count+1"
  26 +done
  27 +
  28 +
  29 +
... ...
batch/RUN_BATCH_Tp.sh 0 → 100755
... ... @@ -0,0 +1,33 @@
  1 +#!/bin/bash
  2 +
  3 +
  4 +name_exp=('tp1' 'tp2' 'tp3' 'tp4' 'tp5'
  5 + 'tp6' 'tp7' 'tp8' 'tp9' 'tp10')
  6 +
  7 +tp="6 7 8 9 10 11 12 13 14 15"
  8 +
  9 +let "count=0"
  10 +for a in $tp
  11 +do
  12 + cat parameter.txt | sed \
  13 + -e "s/nmebatch/${name_exp[$count]}/g" \
  14 + -e "s/probcritbatch/0.15/g" \
  15 + -e "s/strainbatch/3e-5/g" \
  16 + -e "s/Hsbatch/1/g" \
  17 + -e "s/hbatch/4/g" \
  18 + -e "s/Tmbatch/$a/g" \
  19 + > parameter.nml
  20 +
  21 + mv parameter.nml ../nml/
  22 +
  23 +
  24 + qsub ./lanceur.sh
  25 +
  26 + sleep 5s
  27 +
  28 +
  29 + let "count=count+1"
  30 +done
  31 +
  32 +
  33 +
... ...
batch/RUN_BATCH_thickHs.sh 0 → 100755
... ... @@ -0,0 +1,36 @@
  1 +#!/bin/bash
  2 +
  3 +
  4 +h="0.5 1 1.5 2 2.5 3 3.5 4 4.5 5"
  5 +hs="0.2 0.4 0.8 1 1.5 2 2.5 3 3.5 4"
  6 +tp="6 7 8 9 10 11 12 13 14 15"
  7 +let "count=0"
  8 +for a in $h
  9 +do
  10 +for b in $hs
  11 +do
  12 +for c in $tp
  13 +do
  14 + name_exp=thickhstp$count
  15 + cat parameter.txt | sed \
  16 + -e "s/nmebatch/$name_exp/g" \
  17 + -e "s/probcritbatch/0.37/g" \
  18 + -e "s/strainbatch/3e-5/g" \
  19 + -e "s/Hsbatch/$b/g" \
  20 + -e "s/hbatch/$a/g" \
  21 + -e "s/Tmbatch/$c/g" \
  22 + > parameter.nml
  23 +
  24 + mv parameter.nml ../nml/
  25 +
  26 +
  27 + qsub ./lanceur.sh
  28 +
  29 + sleep 30s
  30 +
  31 +
  32 + let "count=count+1"
  33 +done
  34 +done
  35 +done
  36 +
... ...
batch/RUN_BATCH_tpHs.sh 0 → 100755
... ... @@ -0,0 +1,34 @@
  1 +#!/bin/bash
  2 +
  3 +
  4 +tp="6 7 8 9 10 11 12 13 14 15"
  5 +hs="0.2 0.4 0.8 1 1.5 2 2.5 3 3.5 4"
  6 +
  7 +let "count=0"
  8 +for a in $tp
  9 +do
  10 +for b in $hs
  11 +do
  12 + name_exp=exp$count
  13 + cat parameter.txt | sed \
  14 + -e "s/nmebatch/$name_exp/g" \
  15 + -e "s/probcritbatch/0.15/g" \
  16 + -e "s/strainbatch/3e-5/g" \
  17 + -e "s/Hsbatch/$b/g" \
  18 + -e "s/hbatch/4/g" \
  19 + -e "s/Tmbatch/$a/g" \
  20 + > parameter.nml
  21 +
  22 + mv parameter.nml ../nml/
  23 +
  24 +
  25 + qsub ./lanceur.sh
  26 +
  27 + sleep 30s
  28 +
  29 +
  30 + let "count=count+1"
  31 +done
  32 +done
  33 +
  34 +
... ...
batch/lanceur.sh 0 → 100755
... ... @@ -0,0 +1,14 @@
  1 +#!/bin/bash
  2 +#PBS -q default
  3 +#PBS -l walltime=15:00:00
  4 +#PBS -l nodes=1:ppn=20
  5 +#PBS -l nice=19
  6 +#PBS -j oe
  7 +#PBS -N WIM
  8 +module load dot
  9 +hostname
  10 +cd $PBS_O_WORKDIR
  11 +export NPROCS=$(wc -l $PBS_NODEFILE | gawk '//{print $1}')
  12 +cd ..
  13 + ./WIM
  14 +
... ...
batch/parameter.txt 0 → 100644
... ... @@ -0,0 +1,150 @@
  1 +!
  2 +! _______________________
  3 +! | |
  4 +! | WIM PARAMETERS |
  5 +! |_______________________|
  6 +!______________________________________________________________________________
  7 +! WAVES PARAMETERS:
  8 +!
  9 +! Tm -> Peak period [s]
  10 +! Hs -> Significant wave height [m]
  11 +! disp -> Allowing wave dispersion
  12 +! 0: Wave dispersion is not allowed,
  13 +! group speed is the same at all spectrum
  14 +! frequency (cg=max[cg(w)])
  15 +! 1: Wave dispersion is allowed
  16 +!------------------------------------------------------------------------------
  17 +&waves_parameters
  18 +
  19 +Tm =Tmbatch
  20 +Hs =Hsbatch
  21 +disp =0
  22 +
  23 +/
  24 +!______________________________________________________________________________
  25 +! MODEL PARAMETERS:
  26 +!
  27 +! nbin -> Number of grid bin
  28 +! dx -> Spatial resolution [m]
  29 +! Cfl -> Courant–Friedrichs–Lewy condition (0<Cfl<1)
  30 +! Only in the case where disp=0. The CFL condition
  31 +! is needed to calculate the time step.
  32 +! name_sim -> name of the output file
  33 +! root -> destination folder for the output file
  34 +!------------------------------------------------------------------------------
  35 +&model_parameter
  36 +
  37 +nbin =100
  38 +dx =5000
  39 +Cfl =1
  40 +name_sim ='nmebatch'
  41 +root = 'output/'
  42 +
  43 +/
  44 +!______________________________________________________________________________
  45 +! SPECTRUM PARAMETERS:
  46 +!
  47 +! init_spec -> method to build the wave spectrum
  48 +! 2: Swell
  49 +! 1: JONSWAP spectrum
  50 +! 0: Bretschneider spectrum
  51 +! nfreq -> number of frequency bin
  52 +! Tmin -> Minimum period [s]
  53 +! Tmax -> Maximum period [s]
  54 +! alpha_s -> parameter for jonswap spectrum (init_spec=1)
  55 +! beta_s -> parameter for jonswap spectrum (init_spec=1)
  56 +! gamma_s -> parameter for jonswap spectrum (init_spec=1)
  57 +! swell_T -> swell period (init_spec=2)
  58 +! swell_Hs -> swell significant height (init_spec=2)
  59 +!------------------------------------------------------------------------------
  60 +&spectrum_parameters
  61 +
  62 +init_spec =0
  63 +nfreq =800
  64 +Tmin =2.5
  65 +Tmax =20
  66 +alpha_s =0.2044
  67 +beta_s =1.2500
  68 +gamma_s =3.3
  69 +swell_T =19
  70 +swell_Hs=0.09
  71 +
  72 +/
  73 +!______________________________________________________________________________
  74 +! ICE PARAMETERS
  75 +!
  76 +! X_ice -> Distance of the ice edge [m]
  77 +! c_cice -> Ice concentration
  78 +! ice_thick -> method for compute the ice thickness
  79 +! 0: constant thickness
  80 +! 1: thickness is a function of distance
  81 +! from ice edge
  82 +! hice -> Ice thickness (if ice_thick=0) [m]
  83 +! hmax -> Maximum ice thickness (if ice_thick=1) [m]
  84 +! Xh -> distance where h=hmax/2 (if ice_thicl=1) [m]
  85 +! D0 -> initial floe size in the domain [m]
  86 +! gam ->
  87 +! Dmin -> Minimum floe size (if FSD_sheme=1) [m]
  88 +!------------------------------------------------------------------------------
  89 +&ice_parameters
  90 +X_ice =50000
  91 +cice =1
  92 +ice_thick =0
  93 +hice =hbatch
  94 +hmax =hbatch
  95 +Xh =200000
  96 +D0 =500
  97 +gam =1.5
  98 +Dmin =20
  99 +stress_crit =0.67e6
  100 +strain_crit =strainbatch
  101 +P_c =probcritbatch
  102 +
  103 +/
  104 +!________________________________________________________________________________
  105 +! FSD PARAMETERS
  106 +
  107 +! FSD_sheme -> method for compute <D>
  108 +! 0: dumont et al (2011)
  109 +! 1: power law
  110 +!
  111 +! minfloe -> minimum size floe to build the floe size categories [m]
  112 +! maxfloe -> maximum size floe to build the floe size categories [m]
  113 +! nbcat -> number of floe size categories
  114 +!--------------------------------------------------------------------------------
  115 +&fsd_parameters
  116 +FSD_scheme =1
  117 +minfloe =5
  118 +maxfloe =400
  119 +nbcat =60
  120 +
  121 +/
  122 +!_________________________________________________________________________________
  123 +! IDT PARAMETERS
  124 +
  125 +!IDT_scheme -> compute the ice thickness distribution
  126 + 0: no distribution
  127 + 1: distribution (rayleigh)
  128 +
  129 +!mu_IDT -> parameter for the distribution
  130 +!mincat_h -> minimum ice thickness category
  131 +!maxcat_h -> maximum ice thickness category
  132 +!nbcat_h -> number of ice thickness categories
  133 +!---------------------------------------------------------------------------------
  134 +&idt_parameters
  135 +IDT_scheme =1
  136 +mu_IDT =0.5
  137 +mincat_h =0.1
  138 +maxcat_h =10
  139 +nbcat_h =50
  140 +/
  141 +!________________________________________________________________________________
  142 +
  143 +
  144 +
  145 +
  146 +
  147 +
  148 +
  149 +
  150 +
... ...
nml/parameter.nml
... ... @@ -16,7 +16,7 @@
16 16 !------------------------------------------------------------------------------
17 17 &waves_parameters
18 18  
19   -Tm =6
  19 +Tm =8
20 20 Hs =1
21 21 disp =0
22 22  
... ... @@ -31,24 +31,21 @@ disp =0
31 31 ! is needed to calculate the time step.
32 32 ! name_sim -> name of the output file
33 33 ! root -> destination folder for the output file
34   -! FSD_sheme -> method for compute <D>
35   -! 0: dumont et al (2011)
36   -! 1: power law
37 34 !------------------------------------------------------------------------------
38 35 &model_parameter
39 36  
40   -nbin =50
  37 +nbin =100
41 38 dx =5000
42 39 Cfl =1
43   -name_sim ='simulation1'
  40 +name_sim ='temp1'
44 41 root = 'output/'
45   -FSD_scheme =1
46 42  
47 43 /
48 44 !______________________________________________________________________________
49 45 ! SPECTRUM PARAMETERS:
50 46 !
51 47 ! init_spec -> method to build the wave spectrum
  48 +! 2: Swell
52 49 ! 1: JONSWAP spectrum
53 50 ! 0: Bretschneider spectrum
54 51 ! nfreq -> number of frequency bin
... ... @@ -57,16 +54,20 @@ FSD_scheme =1
57 54 ! alpha_s -> parameter for jonswap spectrum (init_spec=1)
58 55 ! beta_s -> parameter for jonswap spectrum (init_spec=1)
59 56 ! gamma_s -> parameter for jonswap spectrum (init_spec=1)
  57 +! swell_T -> swell period (init_spec=2)
  58 +! swell_Hs -> swell significant height (init_spec=2)
60 59 !------------------------------------------------------------------------------
61 60 &spectrum_parameters
62 61  
63   -init_spec =1
64   -nfreq =200
  62 +init_spec =0
  63 +nfreq =800
65 64 Tmin =2.5
66 65 Tmax =20
67 66 alpha_s =0.2044
68 67 beta_s =1.2500
69 68 gamma_s =3.3
  69 +swell_T =19
  70 +swell_Hs=0.09
70 71  
71 72 /
72 73 !______________________________________________________________________________
... ... @@ -87,13 +88,63 @@ gamma_s =3.3
87 88 !------------------------------------------------------------------------------
88 89 &ice_parameters
89 90 X_ice =50000
90   -cice =0.8
91   -ice_thick =1
  91 +cice =1
  92 +ice_thick =0
92 93 hice =2
93   -hmax =4
94   -Xh =100000
  94 +hmax =3
  95 +Xh =200000
95 96 D0 =500
96 97 gam =1.5
97 98 Dmin =20
  99 +stress_crit =0.67e6
  100 +strain_crit =3e-5
  101 +P_c =0.37
  102 +
  103 +/
  104 +!________________________________________________________________________________
  105 +! FSD PARAMETERS
  106 +
  107 +! FSD_sheme -> method for compute <D>
  108 +! 0: dumont et al (2011)
  109 +! 1: power law
  110 +!
  111 +! minfloe -> minimum size floe to build the floe size categories [m]
  112 +! maxfloe -> maximum size floe to build the floe size categories [m]
  113 +! nbcat -> number of floe size categories
  114 +!--------------------------------------------------------------------------------
  115 +&fsd_parameters
  116 +FSD_scheme =1
  117 +minfloe =5
  118 +maxfloe =400
  119 +nbcat =60
  120 +
  121 +/
  122 +!_________________________________________________________________________________
  123 +! IDT PARAMETERS
  124 +
  125 +!IDT_scheme -> compute the ice thickness distribution
  126 + 0: no distribution
  127 + 1: distribution (rayleigh)
  128 +
  129 +!mu_IDT -> parameter for the distribution
  130 +!mincat_h -> minimum ice thickness category
  131 +!maxcat_h -> maximum ice thickness category
  132 +!nbcat_h -> number of ice thickness categories
  133 +!---------------------------------------------------------------------------------
  134 +&idt_parameters
  135 +IDT_scheme =1
  136 +mu_IDT =0.5
  137 +mincat_h =0.1
  138 +maxcat_h =10
  139 +nbcat_h =50
98 140 /
99 141 !________________________________________________________________________________
  142 +
  143 +
  144 +
  145 +
  146 +
  147 +
  148 +
  149 +
  150 +
... ...
output 0 → 120000
... ... @@ -0,0 +1 @@
  1 +/share/work/bauj0001/output
0 2 \ No newline at end of file
... ...
output/data_treatment.m deleted
... ... @@ -1,74 +0,0 @@
1   -clear all;
2   -close all;
3   -
4   -
5   -x=ncread('simulation1.nc','x_axis');
6   -t=ncread('simulation1.nc','time');
7   -om=ncread('simulation1.nc','omega');
8   -spectre=ncread('simulation1.nc','Spectrum');
9   -Dave=ncread('simulation1.nc','Dave');
10   -Dmax=ncread('simulation1.nc','Dmax');
11   -
12   -
13   -f=om/(2*pi);
14   -E=reshape(spectre(end,end,:),length(om),1);
15   -Ei=reshape(spectre(30,:,:),length(x),length(om));
16   -E1=reshape(spectre(1,1,:),length(om),1);
17   -
18   -
19   -
20   -
21   -
22   -figure(1)
23   -cmap=rand(length(t),3);
24   - w = waitforbuttonpress;
25   -for i=1:length(t)
26   -
27   - figure(1)
28   -
29   - sp=reshape(spectre(i,:,:),length(x),length(om));
30   - h=mesh(f,x,sp);
31   - axis([min(f) max(f) min(x) max(x) 0 max(E1)])
32   - xlabel('Frequency [s^{-1}]')
33   - zlabel('Energy')
34   -
35   - ylabel('x [km]')
36   - pause(0.1)
37   -
38   -end
39   -
40   -for i=1:length(om)
41   - EE(i)=sum(Ei(:,i));
42   -
43   -end
44   -figure
45   -plot(E1)
46   -hold on
47   -plot(EE,'r')
48   -
49   -for i=1:length(t)
50   -
51   -E2=reshape(spectre(i,:,:),length(x),length(om));
52   -m0(i)=sum(sum(E2))/sum(E1);
53   -end
54   -figure
55   -plot(t,m0)
56   -
57   -
58   -
59   -
60   -
61   -
62   -figure
63   -
64   -
65   -
66   - plot(x,Dmax,'color','r')
67   - hold on
68   - plot(x,Dave,'--b')
69   - grid on
70   -
71   - xlabel('x [km]')
72   - ylabel('Floe size [m]')
73   -
74   - legend('Dmax','<D>')
75 0 \ No newline at end of file
scripts_matlab/adimentionnal_number.m 0 → 100644
... ... @@ -0,0 +1,69 @@
  1 +
  2 +clear all;
  3 +close all;
  4 +
  5 +
  6 +fig=figure;
  7 +
  8 +% fig.Units='inches';
  9 +% fig.Position=[6 6 9 3.2];
  10 +% fig.OuterPosition=[6 6 9 3.2];
  11 +%
  12 +% fig.PaperUnits = 'inches';
  13 +% fig.PaperPosition = [0 0 9 3.2];
  14 +% fig.PaperPositionMode = 'manual';
  15 +% fig.PaperOrientation='landscape';
  16 +% fig.PaperSize=[9 3.2];
  17 +%
  18 +
  19 +
  20 +for n=3
  21 +titre1=['10^{-6}';'10^{-4}';'10^{-3}';'10^{-2}'];
  22 +titre2=cellstr(titre1);
  23 +
  24 +dos1=['1e6 ';'1e4 ';'test ';'1e2 '];
  25 +dos2=cellstr(dos1);
  26 +
  27 +
  28 +dossier=[num2str(n),'mu1';num2str(n),'mu2';num2str(n),'mu3';dos2(n);num2str(n),'mu4'];
  29 +dossier=cellstr(dossier);
  30 +
  31 +
  32 +diffus=['1e-2_uniforme'];
  33 +
  34 +mix=[1e-6 1e-4 1e-3 1e-2];
  35 +mu=[0.5 1 1.5 2 2.5]/86400;
  36 +
  37 +col=['mo';'ro';'bo';'ko'];
  38 +col=cellstr(col);
  39 +
  40 +for i=1:5
  41 +dos=char(dossier(i));
  42 +
  43 +number(:,n)=sqrt(mu/mix(n));
  44 +
  45 +
  46 +lagbiomass=load(['~/projets/correction_memo/output/',dos,'/weight_',diffus,'.dat']);
  47 +eulbiomass=load(['~/projets/correction_memo/output/',dos,'/eulerian_',diffus,'.dat']);
  48 +
  49 +biolag_tot=sum(lagbiomass,2)*0.5;
  50 +bioeul_tot=sum(eulbiomass,2)*0.5;
  51 +
  52 +std=sqrt((biolag_tot-bioeul_tot).^2./bioeul_tot.^2)*100;
  53 +std(1:3)=0;
  54 +
  55 +maxstd(i,n)=max(std);
  56 +maxbio(i,n)=max(biolag_tot);
  57 +
  58 +
  59 +plot(std)
  60 +hold on
  61 +
  62 +end
  63 +
  64 +end
  65 +
  66 +
  67 +
  68 +
  69 +
... ...
scripts_matlab/anim_joint_distribution.m 0 → 100644
... ... @@ -0,0 +1,99 @@
  1 +clear all;
  2 +close all;
  3 +
  4 +
  5 +dossier=['temp1'];
  6 +path=['/home/bauj0001/projets/WIM/output/',dossier,'/',dossier,'.nc'];
  7 +
  8 +x=ncread(path,'x_axis');
  9 +t=ncread(path,'time');
  10 +om=ncread(path,'omega');
  11 +spectre=ncread(path,'Spectrum');
  12 +Dave=ncread(path,'Dave');
  13 +Dmax=ncread(path,'Dmax');
  14 +FSD=ncread(path,'Floe size distribution');
  15 +Fsize=ncread(path,'floe size');
  16 +thick=ncread(path,'Ice thickness');
  17 +conc=ncread(path,'Ice concentration');
  18 +Hs=ncread(path,'significant height');
  19 +IDT=ncread(path,'Ice Thickness Distribution');
  20 +hcat=ncread(path,'thickness categories');
  21 +
  22 +
  23 +fig=figure(1);
  24 +u = fig.Units;
  25 +fig.Units = 'inches';
  26 +set(fig,'outerposition',[3 3 12 10])
  27 +
  28 +cm=colormap('jet');
  29 +
  30 +
  31 +cm1(1,:)=[1 1 1];
  32 +for i=2:10
  33 + cm1(i,:)=[1-0.08*i 1-0.06*i 1];
  34 +end
  35 +cmap(1:10,:)=cm1;
  36 +cmap(11:59,:)=cm(16:end,:);
  37 +
  38 +
  39 + w = waitforbuttonpress;
  40 +%figure
  41 +%filename = 'joint_distribution_anim.gif';
  42 +
  43 +for j=1:length(t)
  44 +
  45 +FSD1=reshape(FSD(j,j,:,:),length(Fsize),length(hcat));
  46 +
  47 +IDT1=IDT(j,:);
  48 +for i=1:length(hcat)
  49 + FSD1(:,i)=FSD1(:,i)*IDT1(i);
  50 + nfloe(:,i)=((FSD1(:,i)*(5000)^2)./Fsize.^2) ;
  51 +end
  52 + nfloe=nfloe/sum(sum(nfloe));
  53 +
  54 +FSD2(j,:)=sum(FSD1,2);
  55 +subplot(2,2,1)
  56 +hh=surf(hcat,Fsize,nfloe);
  57 + shading interp
  58 +%
  59 +zoom(1.5)
  60 +caxis([0 0.06])
  61 +
  62 +colormap(cmap)
  63 + axis([0.1 10 0 400 0 0.1])
  64 + view(108,60)
  65 +
  66 + xlabel('Ice thickness [m]')
  67 + ylabel('Floe size [m]')
  68 + zlabel('area fraction')
  69 + subplot(2,2,2)
  70 +nfloe2=((FSD2(j,:)*5000^2)./Fsize.^2')/sum((FSD2(j,:)*5000^2)./Fsize.^2');
  71 + plot((Fsize),nfloe2)
  72 +
  73 + axis([0 400 0 1])
  74 + xlabel('Floe size [m]')
  75 + ylabel('normalized number of floe')
  76 +
  77 + subplot(2,2,3:4)
  78 + plot(x,thick)
  79 + hold on
  80 + plot(x(j),thick(j),'ro','markerfacecolor','r')
  81 + hold off
  82 + xlabel('distance [km]')
  83 + ylabel('Average ice thickness [m]')
  84 +%pause(0.02)
  85 + w = waitforbuttonpress;
  86 +
  87 +% drawnow
  88 +% frame = getframe(1);
  89 +% im = frame2im(frame);
  90 +% [A,map] = rgb2ind(im,256);
  91 +% if j == 1;
  92 +% imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',0);
  93 +% else
  94 +% imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0);
  95 +% end
  96 +end
  97 +
  98 +
  99 +
... ...
scripts_matlab/data_treatment.m 0 → 100644
... ... @@ -0,0 +1,154 @@
  1 +clear all;
  2 +close all;
  3 +
  4 +
  5 +x=ncread('simulation1.nc','x_axis');
  6 +t=ncread('simulation1.nc','time');
  7 +om=ncread('simulation1.nc','omega');
  8 +spectre=ncread('simulation1.nc','Spectrum');
  9 +Dave=ncread('simulation1.nc','Dave');
  10 +Dmax=ncread('simulation1.nc','Dmax');
  11 +FSD=ncread('simulation1.nc','Floe size distribution');
  12 +Fsize=ncread('simulation1.nc','floe size');
  13 +thick=ncread('simulation1.nc','Ice thickness');
  14 +conc=ncread('simulation1.nc','Ice concentration');
  15 +Hs=ncread('simulation1.nc','significant height');
  16 +IDT=ncread('simulation1.nc','Ice Thickness Distribution');
  17 +hcat=ncread('simulation1.nc','thickness categories');
  18 +
  19 +
  20 +
  21 +f=om/(2*pi);
  22 +E=reshape(spectre(end,end,:),length(om),1);
  23 +Ei=reshape(spectre(30,:,:),length(x),length(om));
  24 +E1=reshape(spectre(1,1,:),length(om),1);
  25 +
  26 +for i=1:length(f)
  27 + thick2D(:,i)=thick;
  28 +end
  29 +thick2D(:,1)=0;
  30 +thick2D(:,end)=0;
  31 +thick2D(end,:)=0;
  32 +
  33 +
  34 +
  35 +fig=figure(1);
  36 +u = fig.Units;
  37 +fig.Units = 'inches';
  38 +set(fig,'outerposition',[3 3 12 7])
  39 +
  40 +
  41 +cm=colormap('jet');
  42 +
  43 +
  44 +cm1(1,:)=[1 1 1];
  45 +for i=2:10
  46 + cm1(i,:)=[1-0.08*i 1-0.06*i 1];
  47 +end
  48 +cmap(1:10,:)=cm1;
  49 +cmap(11:59,:)=cm(16:end,:);
  50 +
  51 +
  52 +
  53 +filename = 'spectrum_attenuation_anim.gif';
  54 + % w = waitforbuttonpress;
  55 +for i=1:length(t)
  56 +
  57 + figure(1)
  58 +
  59 +
  60 +
  61 + h2=surf(f,x,thick2D*0.05);
  62 + view(-56,22)
  63 +hold on
  64 +
  65 + sp=reshape(spectre(i,:,:),length(x),length(om));
  66 + h=surf(f,x,sp);
  67 + view(-56,22)
  68 + colormap(jet)
  69 + % caxis([-0.1 0.2])
  70 +
  71 + shading interp
  72 + set(h2,'Facealpha',0.3,'facecolor',[0.9 0.9 0.9])
  73 + ll=light;
  74 + set(ll,'position',[0.25 -20 0.4])
  75 + lighting phong
  76 + axis([min(f) max(f) min(x) max(x) 0 max(E1)])
  77 + xlabel('Frequency [s^{-1}]')
  78 + zlabel('Energy')
  79 + hold off
  80 + ylabel('x [km]')
  81 +% pause(0.001)
  82 +%w = waitforbuttonpress;
  83 +drawnow
  84 +frame = getframe(1);
  85 +im = frame2im(frame);
  86 +[A,map] = rgb2ind(im,256);
  87 + if i == 1;
  88 + imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',0);
  89 + else
  90 + imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0);
  91 + end
  92 +
  93 +end
  94 +%
  95 +%
  96 +% for i=1:size(Hs,1)
  97 +%
  98 +% hs2(i)=Hs(i,i);
  99 +% end
  100 +%
  101 +% figure
  102 +% subplot(4,4,1:4)
  103 +% [ax,p1,p2] = plotyy(x,thick,x,conc,'plot','plot');
  104 +%
  105 +% ylabel(ax(1),'Ice thickness [m]')
  106 +% ylabel(ax(2),'Ice concentration')
  107 +% hold on
  108 +% line1=line([55 55], ylim);
  109 +% line2=line([150 150], ylim);
  110 +% line3=line([250 250], ylim);
  111 +% line4=line([400 400], ylim);
  112 +% set(line1,'color','k','linestyle','--')
  113 +% set(line2,'color','k','linestyle','--')
  114 +% set(line3,'color','k','linestyle','--')
  115 +% set(line4,'color','k','linestyle','--')
  116 +%
  117 +%
  118 +% subplot(4,4,5:8)
  119 +%
  120 +% plot(x,hs2)
  121 +%
  122 +% ylabel('Hs [m]')
  123 +%
  124 +%
  125 +% subplot(4,4,9)
  126 +% fsd1=reshape(FSD(:,find(x==55),:),size(FSD,1),length(Fsize));
  127 +% fsd1=fsd1';
  128 +% plot(Fsize,fsd1(:,find(x==55)),'-')
  129 +% xlabel('Floe size [m]')
  130 +% ylabel('area fraction')
  131 +% subplot(4,4,10)
  132 +% fsd1=reshape(FSD(:,find(x==150),:),size(FSD,1),length(Fsize));
  133 +% fsd1=fsd1';
  134 +% plot(Fsize,fsd1(:,find(x==150)),'-')
  135 +% xlabel('Floe size [m]')
  136 +%
  137 +% subplot(4,4,11)
  138 +% fsd1=reshape(FSD(:,find(x==250),:),size(FSD,1),length(Fsize));
  139 +% fsd1=fsd1';
  140 +% plot(Fsize,fsd1(:,find(x==250)),'-')
  141 +% xlabel('Floe size [m]')
  142 +%
  143 +% subplot(4,4,12)
  144 +% fsd1=reshape(FSD(:,find(x==400),:),size(FSD,1),length(Fsize));
  145 +% fsd1=fsd1';
  146 +% plot(Fsize,fsd1(:,find(x==400)),'-')
  147 +% xlabel('Floe size [m]')
  148 +%
  149 +% subplot(4,4,13:16)
  150 +% plot(x,Dave)
  151 +% xlabel('Distance [km]')
  152 +% ylabel('<D> [m]')
  153 +% saveas(gcf,'new_model.png')
  154 +%
... ...
scripts_matlab/fig2_article.m 0 → 100644
... ... @@ -0,0 +1,41 @@
  1 +clear all;
  2 +close all;
  3 +
  4 +
  5 +
  6 +fig=figure;
  7 +
  8 +fig.Units='inches';
  9 +fig.Position=[7.0312 6.6875 4.0312 3.3021];
  10 +fig.OuterPosition=[7.0312 6.6875 4.0312 4.3021];
  11 +
  12 +fig.PaperUnits = 'inches';
  13 +fig.PaperPosition = [0 0 4.0312 3.3021];
  14 +fig.PaperPositionMode = 'manual';
  15 +fig.PaperOrientation='landscape';
  16 +fig.PaperSize=[4.0312 3.3021];
  17 +
  18 +
  19 +strain=[0:1e-6:3e-4];
  20 +prob=exp(-2*3e-5^2./strain.^2);
  21 +
  22 +
  23 + plot(strain(find(prob<0.37)),prob(find(prob<0.37)),'k--','linewidth',1.5)
  24 +hold on
  25 +plot(strain(find(prob>0.37)),prob(find(prob>0.37)),'k-','linewidth',1.5)
  26 +
  27 +a=xlabel('$E_s$');
  28 +b=ylabel('$P_{\varepsilon}=P(E_w>\varepsilon_c)$');
  29 +set(a,'interpreter','latex','fontsize',11)
  30 +set(b,'interpreter','latex','fontsize',11)
  31 +
  32 +scrit=3e-5;
  33 +probcrit=0.37;
  34 +plot(xlim,[probcrit probcrit])
  35 +plot([scrit scrit],ylim)
  36 +c=text(1.5e-4,0.4,'$P_c$');
  37 +d=text(3.5e-5,0.65,'$\varepsilon_c$');
  38 +set(c,'interpreter','latex','fontsize',11)
  39 +set(d,'interpreter','latex','rotation',90,'fontsize',11)
  40 +
  41 +saveas(fig,'Pe_Es','pdf')
... ...
scripts_matlab/fig3_article.m 0 → 100644
... ... @@ -0,0 +1,96 @@
  1 +clear all;
  2 +close all;
  3 +
  4 +
  5 +dossier=['fig1_1'];
  6 +path=['/home/bauj0001/projets/WIM/output/',dossier,'/',dossier,'.nc'];
  7 +
  8 +
  9 +Fsize=ncread(path,'floe size');
  10 +strain1=load('/home/bauj0001/projets/WIM/output/strain1.dat');
  11 +beta1=load('/home/bauj0001/projets/WIM/output/beta1.dat');
  12 +Q1=load('/home/bauj0001/projets/WIM/output/Q1.dat');
  13 +
  14 +strain2=load('/home/bauj0001/projets/WIM/output/strain2.dat');
  15 +beta2=load('/home/bauj0001/projets/WIM/output/beta2.dat');
  16 +Q2=load('/home/bauj0001/projets/WIM/output/Q2.dat');
  17 +
  18 +strain3=load('/home/bauj0001/projets/WIM/output/strain3.dat');
  19 +beta3=load('/home/bauj0001/projets/WIM/output/beta3.dat');
  20 +Q3=load('/home/bauj0001/projets/WIM/output/Q3.dat');
  21 +
  22 +Lmin=pi*0.5*((5e6*2^3)/(3*10*(1-0.3^2)))^0.25;
  23 +Lmin=Lmin*2;
  24 +wavelength=Fsize*2;
  25 +
  26 +
  27 +fig=figure;
  28 +
  29 +fig.Units='inches';
  30 +fig.Position=[7.4167 2.0833 2.9583 5.5312];
  31 +fig.OuterPosition=[7.4167 2.0833 2.9583 6.5312];
  32 +
  33 +fig.PaperUnits = 'inches';
  34 +fig.PaperPosition = [0 0 2.9583 5.5312];
  35 +fig.PaperPositionMode = 'manual';
  36 +fig.PaperOrientation='landscape';
  37 +fig.PaperSize=[2.9583 5.5312];
  38 +
  39 +
  40 +
  41 +
  42 +subplot(3,1,1)
  43 +plot(wavelength(find(wavelength>Lmin)),strain1(find(wavelength>Lmin)),'color',[0.1 0.1 0.1],'linewidth',1.5)
  44 +hold on
  45 +plot(wavelength(find(wavelength>Lmin)),strain2(find(wavelength>Lmin)),'color',[0.5 0.5 0.5],'linewidth',1.5)
  46 +plot(wavelength(find(wavelength>Lmin)),strain3(find(wavelength>Lmin)),'color',[0.8 0.8 0.8],'linewidth',1.5)
  47 +a=ylabel('$E_s$');
  48 +b=xlabel('$\lambda_{\mathrm{ice}}$ [m]');
  49 +set(a,'interpreter','latex')
  50 +set(b,'interpreter','latex')