break_horvat.f90 1.46 KB
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