ice_fracture.f90 5.49 KB
!--------------------------------------------------------------------
!DESCRIPTION: This routine does the wave-induced floe breaking by
!             modifying the floe size distribution.
!--------------------------------------------------------------------


! There is a major issue with as the integral of the distribution
! is not consistent with the ice area and the redistribution is not
! performed conservatively.

subroutine ice_fracture

   use parameters

   !local parameters
   implicit none
        
   !integer                         :: Nband=500
   integer                         :: freq1
   integer                         :: freq2
   integer                         :: ninterp=1

   !double precision                :: TT=13
   double precision                :: m_0
   double precision                :: m_0_stress
   double precision                :: m_0_amp
   double precision                :: m_2_amp
   !double precision                :: Tw
   !double precision                :: kw
   double precision                :: Lmin

   double precision, allocatable   :: P(:)
   double precision, allocatable   :: P_cat(:)
   double precision, allocatable   :: F(:)
   double precision, allocatable   :: Q(:,:)
   double precision, allocatable   :: E_A(:)
   double precision, allocatable   :: Q2(:)
   double precision, allocatable   :: W(:)
   double precision, allocatable   :: strain(:)
   double precision, allocatable   :: stress(:)
   double precision, allocatable   :: wavelength(:)
   double precision, allocatable   :: S_strain(:)
   double precision, allocatable   :: S_stress(:)
   double precision, allocatable   :: k(:)
   double precision, allocatable   :: Dave1(:)
   double precision, allocatable   :: xinterp(:)
   double precision, allocatable   :: Pinterp(:)
   double precision, allocatable   :: nfloe(:,:)
   double precision, allocatable   :: nfloe_tot(:) 
   double precision, allocatable   :: P2(:)

   !allocate(P         (Nband))
   !allocate(S_stress  (Nband))
   !allocate(P2        (Nband))
   allocate(F         (nr   ))
   allocate(P_cat     (nr   ))
   allocate(S_strain  (nr   ))
   allocate(nfloe_tot (nr   ))
   allocate(Q2        (nr   ))
   allocate(Q         (nr,nr))

   allocate(E_A       (nf   ))
   allocate(strain    (nf   ))
   allocate(stress    (nf   ))
   allocate(wavelength(nf   ))
   allocate(k         (nf   ))
   allocate(W         (nf   ))

   allocate(Dave1     (nh   ))
   allocate(nfloe     (nr,nh))

   allocate(xinterp   (ninterp))
   allocate(Pinterp   (ninterp))

   ! calculate the probability of fracture for each parts of the
   ! spectrum and the assiociated wavelength

   do iii=1,nh
      ! Lmin according Mellor 1986
      !Lmin=pi*0.5*((Y*middle_h_cat(iii)**3)/                 &
      !       (3*rhoi*g*(1-nu**2)))**0.25
      ! Lmin corrected from Mellor (1986)
      Lmin=pi*0.25*((Y*middle_h_cat(iii)**3)/                 &
             (3*rhoi*g*(1-nu**2)))**0.25

      W=(g*kice(iii,:))/omega**2
      wavelength=2*pi/kice(iii,:)

      strain=middle_h_cat(iii)*0.5*kice(iii,:)**2*W

      freq1=nf
      freq2=nf

      do ii=1,nr
        
         if ( 2d0*floe_cat(ii+1).lt.wavelength(nf) .or.         &
              2d0*floe_cat(ii+1).gt.wavelength( 1) ) then
            P_cat(ii)=0d0
         else
            do while ( wavelength(freq1).lt.2d0*floe_cat(ii+1) .or. &
                       freq1.lt.0 )       
               freq1=freq1-1
            end do

            ! 0th moment of the strain spectrum
            m_0=sum(E(n,i,freq1:freq2)*strain(freq1:freq2)*domega)
       
            ! significant strain
            S_strain(ii)=2*sqrt(m_0)

            if ( floe_cat(ii).lt.Lmin ) then
               P_cat(ii)=0d0
            else
               P_cat(ii)=exp(-2*strain_crit**2/S_strain(ii)**2)
            end if

            if ( P_cat(ii).lt.P_c ) then
               P_cat(ii)=0d0
            end if 

            freq2=freq1
         end if 
      end do      
     
      P_cat=P_cat*(dx/sum(middle_floe_cat*P_cat+3e-14))
      F(1)=0d0 !the smallest size category cannot break!

      do ii=2,nr
         F(ii)=sum(middle_floe_cat(1:ii)*P_cat(1:ii)/dx)
      end do

      ! this is the redistribution function
      ! the smallest size category cannot redistribute
      Q(:,:)=0
      Q(:,1)=0
      do ii=2,nr
      !Q(1:ii-1,ii)=(middle_floe_cat(2:ii)*P_cat(2:ii))/(sum(middle_floe_cat(2:ii)*P_cat(2:ii)+3e-14))
         Q(1:ii-1,ii)=P_cat(2:ii)/sum(P_cat(2:ii)+3e-14)
      end do

      !compute the new floe size distribution!
      do jj=1,nr
         !Q2(jj)=sum(Q(jj,:)*fsd(n-1,i,:,iii)*F)
         Q2(jj)=sum(Q(jj,:)*fsd(i,:,iii,n-1)*F)
      end do

      !fsd(n,i,1:nr,iii)=fsd(n-1,i,1:nr,iii) -             &
      !                     F*fsd(n-1,i,1:nr,iii) + Q2
      fsd(i,1:nr,iii,n)=fsd(i,1:nr,iii,n-1) -             &
                           F*fsd(i,1:nr,iii,n-1) + Q2
        
      !if (i.eq.25 .and.n.eq.25.and.iii.eq.10)then
      !open(31,file='output/strain3.dat')
      !open(32,file='output/beta3.dat')
      !open(33,file='output/Q3.dat')
        
      !write(31,*)s_strain
      !write(33,*)F
      !do ii=1,nr
      !write(32,*)Q(:,ii)
      !end do
      !end if
   end do

   !compute the new average floe diameter <D>
   do ii=1,nh
      nfloe(:,ii)=fsd(i,:,ii,n)*itd(i,ii)*dx**2/middle_floe_cat**2
      Dave1(ii)=sum(fsd(i,:,ii,n)*middle_floe_cat)
   enddo

   nfloe=nfloe/sum(nfloe+3e-14)

   do ii=1,nr
      nfloe_tot(ii)=sum(nfloe(ii,:))
   end do

   Dave(i)=sum(middle_floe_cat*nfloe_tot)
   !Dave(i)=sum(itd(i,:)*Dave1)

end subroutine ice_fracture