break_horvat.f90 1.46 KB
 Jérémy Baudry committed Nov 23, 2015 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``````