floe_breaking.f90 1.31 KB
subroutine floe_breaking

use parameters

implicit none
	double precision, allocatable :: WN_ice(:),fnt(:),k(:)
	double precision ::rho=1025,d,F,Y=5.5,v=0.3
	integer :: acc=10000
	double precision, allocatable :: W(:)
	double precision, allocatable :: fnt_d(:)
	
	double precision, allocatable :: S(:)
	double precision :: SS
	double precision :: WN_ice_d
	double precision :: RS=7.05e-05
	double precision :: m0_s
	double precision :: m0_a
	double precision :: m2_a
	double precision :: Tw
	double precision :: om_d
	double precision :: wl_w

	allocate(fnt(acc+1))
	allocate(k(acc+1))
	allocate(WN_ice(nfreq))
	allocate(W(nfreq))
	allocate(fnt_d(acc+1))
	allocate(S(nfreq))


	F=Y*h(i)**3/12d0*(1-v**2)
	d=0.9*h(i)

	jj=1
	do ii=0,acc
	k(jj)=real(ii)*10d0/real(acc)
	jj=jj+1
	end do
	
	
	do j=1,nfreq
	fnt=F*k**5+rho*(g-d*omega(j)**2)*k-rho*omega(j)**2
	WN_ice(j)=k(minloc(abs(fnt),DIM=1))
	
	end do

 
	W=g*WN_ice/omega**2
	S=(0.5*h(i))*WN_ice**2*W
	m0_s=sum(E(n,i,1:nfreq)*S**2)*domega
	SS=2*sqrt(m0_s)
	if (SS.gt.RS) then


	m0_a=sum(E(n,i,1:nfreq)*W**2)*domega
	m2_a=sum(omega**2*E(n,i,1:nfreq)*W**2)*domega
	Tw=2*pi*sqrt(m0_a/m2_a)
	om_d=2*pi/Tw

	fnt_d=F*k**5+rho*(g-d*om_d**2)*k-rho*om_d**2
	WN_ice_d=k(minloc(abs(fnt_d),DIM=1))

	wl_w=2*pi/WN_ice_d

	Dmax(i)=max(Dmin,min(0.5*wl_w,Dmax(i)))


	end if




end subroutine floe_breaking