break_horvat.f90 1.46 KB
Newer Older
Jérémy Baudry's avatar
Jérémy Baudry committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
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