Commit c2c1e911c8e0023f290efaa5b3f07ac66f723788

Authored by Dany Dumont
1 parent 963a18cb
Exists in master

Renommage de variables, injection continue de vagues, amelioration du fichier de…

… sortie netcdf, formatage des routines.
batch/parameter.txt
... ... @@ -3,7 +3,7 @@
3 3 ! | |
4 4 ! | WIM PARAMETERS |
5 5 ! |_______________________|
6   -!______________________________________________________________________________
  6 +!____________________________________________________________________
7 7 ! WAVES PARAMETERS:
8 8 !
9 9 ! Tm -> Peak period [s]
... ... @@ -13,7 +13,7 @@
13 13 ! group speed is the same at all spectrum
14 14 ! frequency (cg=max[cg(w)])
15 15 ! 1: Wave dispersion is allowed
16   -!------------------------------------------------------------------------------
  16 +!--------------------------------------------------------------------
17 17 &waves_parameters
18 18  
19 19 Tm =Tmbatch
... ... @@ -21,34 +21,34 @@ Hs =Hsbatch
21 21 disp =0
22 22  
23 23 /
24   -!______________________________________________________________________________
  24 +!____________________________________________________________________
25 25 ! MODEL PARAMETERS:
26 26 !
27   -! nbin -> Number of grid bin
  27 +! nx -> Number of grid bin
28 28 ! dx -> Spatial resolution [m]
29 29 ! Cfl -> Courant–Friedrichs–Lewy condition (0<Cfl<1)
30 30 ! Only in the case where disp=0. The CFL condition
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   -!------------------------------------------------------------------------------
  34 +!--------------------------------------------------------------------
35 35 &model_parameter
36 36  
37   -nbin =100
  37 +nx =100
38 38 dx =5000
39 39 Cfl =1
40 40 name_sim ='nmebatch'
41 41 root = 'output/'
42   -
43 42 /
44   -!______________________________________________________________________________
  43 +
  44 +!___________________________________________________________________
45 45 ! SPECTRUM PARAMETERS:
46 46 !
47 47 ! init_spec -> method to build the wave spectrum
48 48 ! 2: Swell
49 49 ! 1: JONSWAP spectrum
50 50 ! 0: Bretschneider spectrum
51   -! nfreq -> number of frequency bin
  51 +! nf -> number of frequency bin
52 52 ! Tmin -> Minimum period [s]
53 53 ! Tmax -> Maximum period [s]
54 54 ! alpha_s -> parameter for jonswap spectrum (init_spec=1)
... ... @@ -60,7 +60,7 @@ root = &#39;output/&#39;
60 60 &spectrum_parameters
61 61  
62 62 init_spec =0
63   -nfreq =800
  63 +nf =800
64 64 Tmin =2.5
65 65 Tmax =20
66 66 alpha_s =0.2044
... ... @@ -85,8 +85,9 @@ swell_Hs=0.09
85 85 ! D0 -> initial floe size in the domain [m]
86 86 ! gam ->
87 87 ! Dmin -> Minimum floe size (if FSD_sheme=1) [m]
88   -!------------------------------------------------------------------------------
  88 +!--------------------------------------------------------------------
89 89 &ice_parameters
  90 +
90 91 X_ice =50000
91 92 cice =1
92 93 ice_thick =0
... ... @@ -99,7 +100,6 @@ Dmin =20
99 100 stress_crit =0.67e6
100 101 strain_crit =strainbatch
101 102 P_c =probcritbatch
102   -
103 103 /
104 104 !________________________________________________________________________________
105 105 ! FSD PARAMETERS
... ...
nml/parameter.nml
1   -!
2 1 ! _______________________
3 2 ! | |
4 3 ! | WIM PARAMETERS |
5 4 ! |_______________________|
6   -!______________________________________________________________________________
  5 +!____________________________________________________________________
7 6 ! WAVES PARAMETERS:
8 7 !
9 8 ! Tm -> Peak period [s]
... ... @@ -13,42 +12,46 @@
13 12 ! group speed is the same at all spectrum
14 13 ! frequency (cg=max[cg(w)])
15 14 ! 1: Wave dispersion is allowed
16   -!------------------------------------------------------------------------------
  15 +! cont -> Allowing continuous wave input
  16 +! 0: Not continuous (default)
  17 +! 1: Continuous input
  18 +!--------------------------------------------------------------------
17 19 &waves_parameters
18 20  
19   -Tm =8
20   -Hs =1
21   -disp =0
22   -
  21 +Tm =4.5
  22 +Hs =1.5
  23 +disp =1
  24 +cont =1
23 25 /
24   -!______________________________________________________________________________
  26 +
  27 +!____________________________________________________________________
25 28 ! MODEL PARAMETERS:
26 29 !
27   -! nbin -> Number of grid bin
  30 +! nx -> Number of grid bin
28 31 ! dx -> Spatial resolution [m]
29   -! Cfl -> Courant–Friedrichs–Lewy condition (0<Cfl<1)
  32 +! Cfl -> Courant-Friedrich-Lewy condition (0 < Cfl < 1)
30 33 ! Only in the case where disp=0. The CFL condition
31 34 ! is needed to calculate the time step.
32 35 ! name_sim -> name of the output file
33 36 ! root -> destination folder for the output file
34   -!------------------------------------------------------------------------------
  37 +!--------------------------------------------------------------------
35 38 &model_parameter
36 39  
37   -nbin =100
38   -dx =5000
39   -Cfl =1
  40 +nx =100
  41 +dx =5
  42 +Cfl =1.0
40 43 name_sim ='test'
41 44 root = 'output/'
42   -
43 45 /
44   -!______________________________________________________________________________
  46 +
  47 +!____________________________________________________________________
45 48 ! SPECTRUM PARAMETERS:
46 49 !
47 50 ! init_spec -> method to build the wave spectrum
48 51 ! 2: Swell
49 52 ! 1: JONSWAP spectrum
50 53 ! 0: Bretschneider spectrum
51   -! nfreq -> number of frequency bin
  54 +! nf -> number of frequency bin
52 55 ! Tmin -> Minimum period [s]
53 56 ! Tmax -> Maximum period [s]
54 57 ! alpha_s -> parameter for jonswap spectrum (init_spec=1)
... ... @@ -56,21 +59,21 @@ root = &#39;output/&#39;
56 59 ! gamma_s -> parameter for jonswap spectrum (init_spec=1)
57 60 ! swell_T -> swell period (init_spec=2)
58 61 ! swell_Hs -> swell significant height (init_spec=2)
59   -!------------------------------------------------------------------------------
  62 +!--------------------------------------------------------------------
60 63 &spectrum_parameters
61 64  
62 65 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   -
  66 +nf =36
  67 +Tmin =2.5
  68 +Tmax =20
  69 +alpha_s =0.2044
  70 +beta_s =1.2500
  71 +gamma_s =3.3
  72 +swell_T =19
  73 +swell_Hs =0.09
72 74 /
73   -!______________________________________________________________________________
  75 +
  76 +!____________________________________________________________________
74 77 ! ICE PARAMETERS
75 78 !
76 79 ! X_ice -> Distance of the ice edge [m]
... ... @@ -82,62 +85,57 @@ swell_Hs=0.09
82 85 ! hice -> Ice thickness (if ice_thick=0) [m]
83 86 ! hmax -> Maximum ice thickness (if ice_thick=1) [m]
84 87 ! Xh -> distance where h=hmax/2 (if ice_thicl=1) [m]
85   -!------------------------------------------------------------------------------
  88 +!--------------------------------------------------------------------
86 89 &ice_parameters
87   -X_ice =50000
88   -cice =1
89   -ice_thick =0
90   -hice =3
91   -hmax =3
92   -Xh =200000
93   -strain_crit =3e-5
94   -P_c =0.55
95 90  
  91 +X_ice =10
  92 +cice =1
  93 +ice_thick=0
  94 +hice =0
  95 +hmax =0.8
  96 +Xh =250
  97 +strain_crit=3e-5
  98 +P_c =0.55
96 99 /
97   -!________________________________________________________________________________
  100 +
  101 +!____________________________________________________________________
98 102 ! FSD PARAMETERS
99 103  
100 104 ! FSD_sheme -> method for compute <D>
101 105 ! 0: dumont et al (2011)
102 106 ! 1: power law
103 107 !
104   -! minfloe -> minimum size floe to build the floe size categories [m]
105   -! maxfloe -> maximum size floe to build the floe size categories [m]
106   -! nbcat -> number of floe size categories
107   -!--------------------------------------------------------------------------------
  108 +! minfloe -> minimum floe size of the FSD [m]
  109 +! maxfloe -> maximum floe size of the FSD [m]
  110 +! nr -> number of floe size categories
  111 +!--------------------------------------------------------------------
108 112 &fsd_parameters
  113 +
109 114 FSD_scheme =1
110 115 minfloe =5
111 116 maxfloe =400
112   -nbcat =60
113   -
  117 +nr =60
114 118 /
115   -!_________________________________________________________________________________
116   -! IDT PARAMETERS
117 119  
118   -!IDT_scheme -> compute the ice thickness distribution
  120 +!____________________________________________________________________
  121 +! ITD PARAMETERS
  122 +
  123 +!itd_scheme -> compute the ice thickness distribution
119 124 0: no distribution
120 125 1: distribution (rayleigh)
121 126  
122   -!mu_IDT -> parameter for the distribution
  127 +!mu_itd -> parameter for the distribution
123 128 !mincat_h -> minimum ice thickness category
124 129 !maxcat_h -> maximum ice thickness category
125   -!nbcat_h -> number of ice thickness categories
126   -!---------------------------------------------------------------------------------
127   -&idt_parameters
128   -IDT_scheme =1
129   -mu_IDT =0.5
  130 +!nh -> number of ice thickness categories
  131 +!--------------------------------------------------------------------
  132 +&itd_parameters
  133 +
  134 +itd_scheme =1
  135 +mu_itd =0.5
130 136 mincat_h =0.1
131 137 maxcat_h =10
132   -nbcat_h =50
  138 +nh =50
133 139 /
134   -!________________________________________________________________________________
135   -
136   -
137   -
138   -
139   -
140   -
141   -
142   -
143 140  
  141 +!____________________________________________________________________
... ...
src/advection.f90
  1 +!-------------------------------------------------------------------
  2 +!DESCRIPTION: Here, the wave spectrum is advected through the
  3 +! domain. The advection equation is solved using a
  4 +! Lax-wendroff discretization sheme with a Superbee
  5 +! flux limiter.
  6 +!-------------------------------------------------------------------
1 7  
2   -!_______________________________________________________________________________
  8 +!INTERFACE:
  9 +subroutine advection
3 10  
4   - !DESCRIPTION: Here, the wave spectrum is advected through the
5   - !domain. The advection equation is solved using a Lax-wendroff
6   - !discretization sheme with a superbee flux limiter.
  11 + !MODULE USES:
  12 + use parameters
7 13  
8   -!_______________________________________________________________________________
9   -
10   -
11   - !INTERFACE:
12   - subroutine advection
13   -
14   - !MODULE USES:
15   - use parameters
16   -
17   - !LOCAL PARAMETERS:
18   - implicit none
  14 + !LOCAL PARAMETERS:
  15 + implicit none
19 16  
20   - double precision, allocatable :: diffF(:)
21   - double precision, allocatable :: diffl(:)
22   - double precision, allocatable :: diffr(:)
23   - double precision, allocatable :: phi(:)
24   - double precision, allocatable :: theta(:)
25   - double precision, allocatable :: F(:)
26   - double precision, allocatable :: diff1(:)
27   -
  17 + double precision, allocatable :: diffF(:)
  18 + double precision, allocatable :: diffl(:)
  19 + double precision, allocatable :: diffr(:)
  20 + double precision, allocatable :: phi (:)
  21 + double precision, allocatable :: theta(:)
  22 + double precision, allocatable :: F (:)
  23 + double precision, allocatable :: diff1(:)
28 24  
29   - allocate(diffl(nbin))
30   - allocate(diffr(nbin))
31   - allocate(diffF(nbin))
32   - allocate(theta(nbin))
33   - allocate(phi(nbin))
34   - allocate(F(nbin))
35   - allocate(diff1(nbin-1))
36   -
37   -!________________________________________________________________________________
38   -
39   - do i=2,nbin
40   - diffl(i)=E(n-1,i,ii)-E(n-1,i-1,ii)
41   - end do
42   -
43   - do i=1,nbin-1
44   - diffr(i)=E(n-1,i+1,ii)-E(n-1,i,ii)
45   - end do
46   -
47   - diffr(nbin)=-E(n-1,nbin,ii) !set dE=0 at boundaries
48   - diffl(1)=E(n-1,1,ii)
49   -
50   -
51   - do i=1,nbin !Superbee flux limiter
52   - theta(i)=diffl(i)/(diffr(i)+3e-14)
53   - phi(i)=max(0d0,min(theta(i)*2,1d0),min(theta(i),2d0))
54   - end do
55   -
56   -
57   - !Lax-Wendroff sheme:
58   - F=E(n-1,1:nbin,ii)+0.5*(1-CN(ii))*diffr*phi
59   - do i=2,nbin
60   - diffF(i)=F(i)-F(i-1)
61   - end do
62   - diffF(1)=F(1)
63   - E(n,1:nbin,ii)=E(n-1,1:nbin,ii)-CN(ii)*diffF
64   -
65   -
66   -!________________________________________________________________________________
  25 + allocate(diffl(nx ))
  26 + allocate(diffr(nx ))
  27 + allocate(diffF(nx ))
  28 + allocate(theta(nx ))
  29 + allocate(phi (nx ))
  30 + allocate(F (nx ))
  31 + allocate(diff1(nx-1))
  32 +
  33 + if ( cont.eq.1 ) then
  34 + E(1,ii,n-1)=Ei(ii)
  35 + end if
  36 +
  37 + do i=2,nx
  38 + diffl(i) = E(i ,ii,n-1) - E(i-1,ii,n-1)
  39 + end do
  40 +
  41 + do i=1,nx-1
  42 + diffr(i) = E(i+1,ii,n-1) - E(i ,ii,n-1)
  43 + end do
  44 +
  45 + ! set dE=0 at boundaries
  46 + diffr(nx)=-E(nx,ii,n-1)
  47 + diffl(1 )= E( 1,ii,n-1)
  48 +
  49 + ! Superbee flux limiter
  50 + do i=1,nx
  51 + theta(i)=diffl(i)/(diffr(i)+3e-14)
  52 + phi(i)=max(0d0,min(theta(i)*2,1d0),min(theta(i),2d0))
  53 + end do
  54 +
  55 + ! Lax-Wendroff sheme
  56 + F=E(1:nx,ii,n-1) + 0.5*(1-CN(ii))*diffr*phi
  57 +
  58 + do i=2,nx
  59 + diffF(i)=F(i) - F(i-1)
  60 + end do
  61 +
  62 + diffF(1)=F(1)
  63 + E(1:nx,ii,n)=E(1:nx,ii,n-1)-CN(ii)*diffF
  64 +
67 65 end subroutine advection
68   -
69   -
70   -
71   -
72   -
73   -
74   -
75   -
76   -
77   -
78   -
79   -
80   -
81   -
82   -
83   -
... ...
src/attenuation.f90
1 1  
2   -!____________________________________________________________________________
  2 +!____________________________________________________________________
  3 +!
  4 +!DESCRIPTION: In this routine, the attenuation coefficient is
  5 +! calculated according Kohout and Meylan (2008) and the
  6 +! spectrum attenuation is computed.
  7 +!
  8 +!____________________________________________________________________
3 9  
  10 +!INTERFACE:
  11 +subroutine attenuation
4 12  
5   - !DESCRIPTION: In this routine, the attenuation
6   - !coefficient is calculated according Kohout and
7   - !Meylan (2008) and the spectrum attenuation is
8   - !computed.
  13 +!USES:
  14 +use parameters
9 15  
10   -!____________________________________________________________________________
  16 +implicit none
11 17  
  18 +double precision :: q11,q12,q13,q14,q15,q21,q22,q23,q24,q25
  19 +double precision :: q31,q32,q33,q34,q35
  20 +double precision, allocatable :: p1(:),p2(:),p3(:)
  21 +double precision, allocatable :: att(:)
12 22  
  23 +allocate(p1 (nf))
  24 +allocate(p2 (nf))
  25 +allocate(p3 (nf))
  26 +allocate(att(nf))
13 27  
14   - !INTERFACE:
15   - subroutine attenuation
16   -
17   - !USES:
18   - use parameters
19   -
20   - implicit none
21   -
22   - double precision :: q11,q12,q13,q14,q15,q21,q22,q23,q24,q25
23   - double precision :: q31,q32,q33,q34,q35
24   - double precision, allocatable :: p1(:),p2(:),p3(:)
25   - double precision, allocatable :: att(:)
26   -
27   - allocate(p1(nfreq))
28   - allocate(p2(nfreq))
29   - allocate(p3(nfreq))
30   - allocate(att(nfreq))
31   -
32   -!____________________________________________________________________________
33 28  
34 29 q11 = -0.00000777
35 30 q12 = 0.00032080
... ... @@ -53,21 +48,17 @@ p1 = q11*T**4 + q12*T**3 + q13*T**2 + q14*T + q15
53 48 p2 = q21*T**4 + q22*T**3 + q23*T**2 + q24*T + q25
54 49 p3 = q31*T**4 + q32*T**3 + q33*T**2 + q34*T + q35
55 50  
56   -alpha(i,1:nfreq)=-1*(p1*h(i)**2 + p2*h(i) + p3)
57   -do j=1,nfreq
58   -if (alpha(i,j).lt.0)then
59   -alpha(i,j)=0d0
60   -end if
61   -end do
62   -
63   -
64   -att=C_ice(i)*alpha(i,1:nfreq)/(Dave(i)+3e-14)
65   -
66   -S_ice(n,i,1:nfreq)=-att*E(n,i,1:nfreq)
67   -
68   -E(n,i,1:nfreq)=E(n,i,1:nfreq)*exp(-att*Cg*dt)
  51 +alpha(i,:)=-1*(p1*h(i)**2 + p2*h(i) + p3)
69 52  
  53 +do ii=1,nf
  54 + if ( alpha(i,ii).lt.0 )then
  55 + alpha(i,ii)=0d0
  56 + end if
  57 +end do
70 58  
71   -!_____________________________________________________________________________
  59 +att=C_ice(i)*alpha(i,:)/(Dave(i)+3e-14)
  60 +!S_ice(n,i,1:nf)=-att*E(n,i,1:nf)
  61 +S_ice(i,:,n)=-att*E(i,:,n)
  62 +E(i,:,n)=E(i,:,n)*exp(-att*Cg*dt)
72 63  
73 64 end subroutine attenuation
... ...
src/ice_fracture.f90
1   - subroutine ice_fracture
2 1  
3   - use parameters
  2 +subroutine ice_fracture
4 3  
5   - !local parameters
6   - implicit none
7   -
8   - double precision :: Y=5.5e9
9   - double precision :: TT=13
10   - double precision :: poisson=0.3
11   - double precision, allocatable :: P(:)
12   - double precision, allocatable :: P_cat(:)
13   - double precision, allocatable :: F(:)
14   - double precision, allocatable :: Q(:,:)
15   - double precision, allocatable :: E_A(:)
16   - double precision, allocatable :: Q2(:)
17   - double precision, allocatable :: W(:)
18   - double precision, allocatable :: strain(:)
19   - double precision, allocatable :: stress(:)
20   - integer :: Nband=500
21   - double precision :: m_0
22   - double precision :: m_0_stress
23   - double precision :: m_0_amp
24   - double precision :: m_2_amp
25   - double precision :: Tw
26   - double precision :: kw
27   - double precision,allocatable :: wavelength(:)
28   - double precision,allocatable :: S_strain(:)
29   - double precision, allocatable :: S_stress(:)
30   - integer :: freq1
31   - integer :: freq2
32   - double precision, allocatable :: k(:)
33   - double precision :: Lmin
34   - double precision, allocatable :: Dave1(:)
35   - integer ::ninterp=1
36   - double precision, allocatable :: xinterp(:)
37   - double precision, allocatable :: Pinterp(:)
38   - double precision, allocatable ::nfloe(:,:)
39   - double precision, allocatable :: nfloe_tot(:)
40   - double precision, allocatable :: P2(:)
41   -
42   - allocate(P(Nband))
43   - allocate(P_cat(nbcat))
44   - allocate(F(nbcat))
45   - allocate(Q(nbcat,nbcat))
46   - allocate(E_A(nfreq))
47   - allocate(Q2(nbcat))
48   - allocate(strain(nfreq))
49   - allocate(stress(nfreq))
50   - allocate(wavelength(nfreq))
51   - allocate(S_strain(nbcat))
52   - allocate(k(nfreq))
53   - allocate(S_stress(Nband))
54   - allocate(Dave1(nbcat_h))
55   - allocate(xinterp(ninterp))
56   - allocate(Pinterp(ninterp))
57   - allocate(nfloe(nbcat,nbcat_h))
58   - allocate(nfloe_tot(nbcat))
59   - allocate(P2(Nband))
60   - allocate(W(nfreq))
61   -
62   -
63   -
64   - !calculate the probability of fracture for each parts of the
65   - !spectrum and the assiociated wavelength
66   -
67   -
68   - do iii=1,nbcat_h
69   -
70   - Lmin=pi*0.5*((5e6*middle_h_cat(iii)**3)/(3*10*(1-0.3**2)))**0.25!Lmin according Mellor 1986
71   -
72   -
73   - W=(9.81*kice(iii,:))/omega**2
74   - wavelength=2*pi/kice(iii,:)
  4 + use parameters
75 5  
76   - strain=middle_h_cat(iii)*0.5*kice(iii,:)**2*W
77   -
78   - freq1=nfreq
79   - freq2=nfreq
80   -
81   - do ii=1,nbcat
  6 + !local parameters
  7 + implicit none
82 8  
83   - if(2d0*floe_cat(ii+1).lt.wavelength(nfreq).or.2d0*floe_cat(ii+1).gt.wavelength(1)) then
84   - P_cat(ii)=0d0
85   - else
86   -
87   -
88   - do while(wavelength(freq1).lt.2d0*floe_cat(ii+1).or.freq1.lt.0)
89   - freq1=freq1-1
90   - end do
91   -
92   - !0th moment of the strain spectrum
93   - m_0=sum(E(n,i,freq1:freq2)*strain(freq1:freq2)*domega)
94   -
95   - !significant strain
96   - S_strain(ii)=2*sqrt(m_0)
  9 + integer :: Nband=500
  10 + integer :: freq1
  11 + integer :: freq2
  12 + integer :: ninterp=1
  13 +
  14 + double precision :: Y=5.5e9
  15 + double precision :: TT=13
  16 + double precision :: poisson=0.3
  17 + double precision :: m_0
  18 + double precision :: m_0_stress
  19 + double precision :: m_0_amp
  20 + double precision :: m_2_amp
  21 + double precision :: Tw
  22 + double precision :: kw
  23 + double precision :: Lmin
  24 +
  25 + double precision, allocatable :: P(:)
  26 + double precision, allocatable :: P_cat(:)
  27 + double precision, allocatable :: F(:)
  28 + double precision, allocatable :: Q(:,:)
  29 + double precision, allocatable :: E_A(:)
  30 + double precision, allocatable :: Q2(:)
  31 + double precision, allocatable :: W(:)
  32 + double precision, allocatable :: strain(:)
  33 + double precision, allocatable :: stress(:)
  34 + double precision, allocatable :: wavelength(:)
  35 + double precision, allocatable :: S_strain(:)
  36 + double precision, allocatable :: S_stress(:)
  37 + double precision, allocatable :: k(:)
  38 + double precision, allocatable :: Dave1(:)
  39 + double precision, allocatable :: xinterp(:)
  40 + double precision, allocatable :: Pinterp(:)
  41 + double precision, allocatable :: nfloe(:,:)
  42 + double precision, allocatable :: nfloe_tot(:)
  43 + double precision, allocatable :: P2(:)
  44 +
  45 + !allocate(P (Nband))
  46 + !allocate(S_stress (Nband))
  47 + !allocate(P2 (Nband))
  48 + allocate(F (nr ))
  49 + allocate(P_cat (nr ))
  50 + allocate(S_strain (nr ))
  51 + allocate(nfloe_tot (nr ))
  52 + allocate(Q2 (nr ))
  53 + allocate(Q (nr,nr))
  54 +
  55 + allocate(E_A (nf ))
  56 + allocate(strain (nf ))
  57 + allocate(stress (nf ))
  58 + allocate(wavelength(nf ))
  59 + allocate(k (nf ))
  60 + allocate(W (nf ))
  61 +
  62 + allocate(Dave1 (nh ))
  63 + allocate(nfloe (nr,nh))
  64 +
  65 + allocate(xinterp (ninterp))
  66 + allocate(Pinterp (ninterp))
  67 +
  68 + ! calculate the probability of fracture for each parts of the
  69 + ! spectrum and the assiociated wavelength
  70 +
  71 + do iii=1,nh
  72 + ! Lmin according Mellor 1986
  73 + Lmin=pi*0.5*((5e6*middle_h_cat(iii)**3)/(3*10*(1-0.3**2)))**0.25
97 74  
98   -
99   -
100   - if(floe_cat(ii).lt.Lmin)then
101   - P_cat(ii)=0d0
102   - else
103   - P_cat(ii)=exp(-2*strain_crit**2/S_strain(ii)**2)
104   - end if
  75 + W=(9.81*kice(iii,:))/omega**2
  76 + wavelength=2*pi/kice(iii,:)
105 77  
106   - if (P_cat(ii).lt.P_c) then
107   - P_cat(ii)=0d0
108   - end if
  78 + strain=middle_h_cat(iii)*0.5*kice(iii,:)**2*W
109 79  
  80 + freq1=nf
  81 + freq2=nf
110 82  
111   -
112   - freq2=freq1
113   - end if
114   - end do
  83 + do ii=1,nr
115 84  
116   -
  85 + if ( 2d0*floe_cat(ii+1).lt.wavelength(nf) .or. &
  86 + 2d0*floe_cat(ii+1).gt.wavelength(1)) then
  87 + P_cat(ii)=0d0
  88 + else
  89 + do while ( wavelength(freq1).lt.2d0*floe_cat(ii+1) .or. &
  90 + freq1.lt.0 )
  91 + freq1=freq1-1
  92 + end do
  93 +
  94 + ! 0th moment of the strain spectrum
  95 + m_0=sum(E(n,i,freq1:freq2)*strain(freq1:freq2)*domega)
  96 +
  97 + ! significant strain
  98 + S_strain(ii)=2*sqrt(m_0)
  99 +
  100 + if ( floe_cat(ii).lt.Lmin ) then
  101 + P_cat(ii)=0d0
  102 + else
  103 + P_cat(ii)=exp(-2*strain_crit**2/S_strain(ii)**2)
  104 + end if
  105 +
  106 + if ( P_cat(ii).lt.P_c ) then
  107 + P_cat(ii)=0d0
  108 + end if
  109 +
  110 + freq2=freq1
  111 + end if
  112 + end do
117 113  
118   - P_cat=P_cat*(dx/sum(middle_floe_cat*P_cat+3e-14))
119   - F(1)=0d0 !the smallest size category cannot break!
120   -
121   - do ii=2,nbcat
122   - F(ii)=sum(middle_floe_cat(1:ii)*P_cat(1:ii)/dx)
123   - end do
124   -
125   -
126   -
127   - Q(:,:)=0
128   - !this is the redistribution function
129   - Q(:,1)=0 ! the smallest size category cannot redistribute
130   - do ii=2,nbcat
131   - !Q(1:ii-1,ii)=(middle_floe_cat(2:ii)*P_cat(2:ii))/(sum(middle_floe_cat(2:ii)*P_cat(2:ii)+3e-14))
132   - Q(1:ii-1,ii)=P_cat(2:ii)/sum(P_cat(2:ii)+3e-14)
133   - end do
134   -
135   -
136   -
137   - !compute the new floe size distribution!
138   - do jj=1,nbcat
139   -
140   - Q2(jj)=sum(Q(jj,:)*FSD(n-1,i,:,iii)*F)
141   - end do
142   -
143   - FSD(n,i,1:nbcat,iii)=FSD(n-1,i,1:nbcat,iii)-F*FSD(n-1,i,1:nbcat,iii)+Q2
144   -
  114 + P_cat=P_cat*(dx/sum(middle_floe_cat*P_cat+3e-14))
  115 + F(1)=0d0 !the smallest size category cannot break!
  116 +
  117 + do ii=2,nr
  118 + F(ii)=sum(middle_floe_cat(1:ii)*P_cat(1:ii)/dx)
  119 + end do
  120 +
  121 + ! this is the redistribution function
  122 + ! the smallest size category cannot redistribute
  123 + Q(:,:)=0
  124 + Q(:,1)=0
  125 + do ii=2,nr
  126 + !Q(1:ii-1,ii)=(middle_floe_cat(2:ii)*P_cat(2:ii))/(sum(middle_floe_cat(2:ii)*P_cat(2:ii)+3e-14))
  127 + Q(1:ii-1,ii)=P_cat(2:ii)/sum(P_cat(2:ii)+3e-14)
  128 + end do
  129 +
  130 + !compute the new floe size distribution!
  131 + do jj=1,nr
  132 + Q2(jj)=sum(Q(jj,:)*FSD(n-1,i,:,iii)*F)
  133 + end do
  134 +
  135 + FSD(n,i,1:nr,iii)=FSD(n-1,i,1:nr,iii) - &
  136 + F*FSD(n-1,i,1:nr,iii) + Q2
145 137  
146   - !if (i.eq.25 .and.n.eq.25.and.iii.eq.10)then
147   - !open(31,file='output/strain3.dat')
148   - !open(32,file='output/beta3.dat')
149   - !open(33,file='output/Q3.dat')
  138 + !if (i.eq.25 .and.n.eq.25.and.iii.eq.10)then
  139 + !open(31,file='output/strain3.dat')
  140 + !open(32,file='output/beta3.dat')
  141 + !open(33,file='output/Q3.dat')
150 142  
151   - !write(31,*)s_strain
152   - !write(33,*)F
153   - !do ii=1,nbcat
154   - !write(32,*)Q(:,ii)
155   - !end do
156   - !end if
157   - end do
158   - !compute the new average floe diameter <D>
159   - do ii=1,nbcat_h
160   - nfloe(:,ii)=FSD(n,i,:,ii)*IDT(i,ii)*dx**2/middle_floe_cat**2
161   - Dave1(ii)=sum(FSD(n,i,:,ii)*middle_floe_cat)
162   - enddo
163   - nfloe=nfloe/sum(nfloe+3e-14)
164   - do ii=1,nbcat
165   - nfloe_tot(ii)=sum(nfloe(ii,:))
166   - end do
167   - Dave(i)=sum(middle_floe_cat*nfloe_tot)
168   - !Dave(i)=sum(IDT(i,:)*Dave1)
169   - end subroutine
170   -
171   -
172   -
  143 + !write(31,*)s_strain
  144 + !write(33,*)F
  145 + !do ii=1,nr
  146 + !write(32,*)Q(:,ii)
  147 + !end do
  148 + !end if
  149 + end do
  150 +
  151 + !compute the new average floe diameter <D>
  152 + do ii=1,nh
  153 + nfloe(:,ii)=FSD(n,i,:,ii)*itd(i,ii)*dx**2/middle_floe_cat**2
  154 + Dave1(ii)=sum(FSD(n,i,:,ii)*middle_floe_cat)
  155 + enddo
  156 +
  157 + nfloe=nfloe/sum(nfloe+3e-14)
  158 +
  159 + do ii=1,nr
  160 + nfloe_tot(ii)=sum(nfloe(ii,:))
  161 + end do
  162 +
  163 + Dave(i)=sum(middle_floe_cat*nfloe_tot)
  164 + !Dave(i)=sum(itd(i,:)*Dave1)
  165 +
  166 +end subroutine ice_fracture
... ...
src/initialization.f90
1 1 !
2   -!_______________________________________________________________________
  2 +!____________________________________________________________________
3 3 !
4   -! DESCRIPTION:
5   -! This is the initialization routine which contains the
6   -! initial spectrum construction, set initial values for arrays
7   -! and construct the ice transect.
8   -!_______________________________________________________________________
9   -
10   - !INTERFACE:
11   - subroutine initialization
12   -
13   - !MODULE USES:
14   - use parameters
15   -
16   - !LOCAL PARAMETERS AND VARIABLES
17   - implicit none
18   - double precision, allocatable ::Gf(:),PM(:)
19   - integer ::X1
20   - complex(kind=8), dimension(6) :: poly
21   - complex(kind=8), dimension(5) :: roots
22   - logical :: param1
23   - logical :: param2
24   - allocate(Gf(nfreq))
25   - allocate(PM(nfreq))
26   -
27   -
28   -
29   -
30   -
31   - !construct time array
32   - time(1)=0
33   - do ii=2,nsteps
34   - time(ii)=time(ii-1)+dt/60
35   - end do
36   -
37   - !construct space array
38   - x_axis(1)=0
39   - do ii=2,nbin
40   - x_axis(ii)=x_axis(ii-1)+dx/1000
41   - end do
42   -
43   -
44   -
45   -!_________________________INITIAL SPECTRUM_____________________________
46   -
47   - E(1:nsteps,1:nbin,1:nfreq)=0d0 !INITIALIZE SPECTRUM ARRAY
48   -
49   - if(init_spec.eq.1) then !use JONSWAP spectrum
50   -
51   - do i=1,nfreq
52   - if (freq(i).le.freq_s) then
53   - sigma_s(i)=0.07
54   - else
55   - sigma_s(i)=0.09
56   - end if
57   - end do
58   -
59   - Gf=gamma_s**(exp((-(freq-freq_s)**2)/(2*sigma_s**2*freq_s**2)))
60   - PM=alpha_s*Hs**2*(freq_s**4/freq**5)*exp(-beta_s*(freq_s/freq)**4)
61   -
62   - Ei=Gf*PM
  4 +!DESCRIPTION: This is the initialization routine which contains the
  5 +! initial spectrum construction, set initial values for
  6 +! arrays and construct the ice transect.
  7 +!____________________________________________________________________
  8 +
  9 +!INTERFACE:
  10 +subroutine initialization
  11 +
  12 +!MODULE USES:
  13 + use parameters
  14 +
  15 +!LOCAL PARAMETERS AND VARIABLES
  16 + implicit none
  17 +
  18 + double precision, allocatable :: Gf(:),PM(:)
  19 +
  20 + integer :: x1
  21 +
  22 + complex(kind=8), dimension(6) :: poly
  23 + complex(kind=8), dimension(5) :: roots
  24 +
  25 + logical :: param1
  26 + logical :: param2
  27 +
  28 + allocate(Gf(nf))
  29 + allocate(PM(nf))
  30 +
  31 +! construct time array in hours
  32 + time(1)=0
  33 + do ii=2,nt
  34 + time(ii)=time(ii-1)+dt/3600
  35 + end do
  36 +
  37 + write(*,*) ' * Time step is ',dt,' s'
  38 +
  39 +! construct space array in km
  40 + x_axis(1)=0
  41 + do ii=2,nx
  42 + x_axis(ii)=x_axis(ii-1)+dx/1000
  43 + end do
  44 +
  45 +
  46 +!--------------------- INITIAL SPECTRUM -----------------------------
  47 +
  48 + E(:,:,:)=0d0
  49 +
  50 + if ( init_spec.eq.1 ) then
  51 +
  52 + ! JONSWAP spectrum
  53 + do i=1,nf
  54 + if ( freq(i).le.freq_s ) then
  55 + sigma_s(i)=0.07
  56 + else
  57 + sigma_s(i)=0.09
  58 + end if
  59 + end do
  60 +
  61 + Gf=gamma_s**(exp((-(freq-freq_s)**2)/(2*sigma_s**2*freq_s**2)))
  62 + PM=alpha_s*Hs**2*(freq_s**4/freq**5)*exp(-beta_s*(freq_s/freq)**4)
  63 + Ei=Gf*PM
63 64  
64   - else if (init_spec.eq.0) then !use bretschneider spectrum
65   -
66   - Ei=(1.25*Hs**2*(1/freq)**5)/(8*pi*Tm**4)*exp(-1.25*((1/freq)/Tm)**4)
67   -
68   - else
69   - Ei=1/(0.01*sqrt(2*pi))*exp(-(omega-2*pi/swell_T)**2/(2*0.01**2))
70   - Ei=(swell_Hs/4d0)**2*Ei/(sum(Ei*domega))
71   - end if
72   -
  65 + else if ( init_spec.eq.0 ) then
73 66  
  67 + ! Bretschneider spectrum
  68 + Ei=(1.25*Hs**2*(1/freq)**5)/(8*pi*Tm**4)* &
  69 + exp(-1.25*((1/freq)/Tm)**4)
74 70  
75   - E(1,1,1:nfreq)=Ei !initial spectrum
  71 + else
76 72  
77   -!_______________________________________________________________________
  73 + Ei=1/(0.01*sqrt(2*pi))*exp(-(omega-2*pi/swell_T)**2/(2*0.01**2))
  74 + Ei=(swell_Hs/4d0)**2*Ei/(sum(Ei*domega))
78 75  
79   -!_______________________ICE_TRANSECT____________________________________
80   -S_ice(1,1,1:nfreq)=0
81   -X1=floor(X_ice/dx) !find in which grid bin is the ice edge
82   -C_ice(1:X1)=0 !ice concentration is 0 before ice edge!
83   -C_ice(X1:nbin)=cice !ice concentration in the transect