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