ice_fracture.f90 4.84 KB

subroutine ice_fracture

   use parameters

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

   double precision                :: Y=5.5e9
   double precision                :: TT=13
   double precision                :: poisson=0.3
   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*((5e6*middle_h_cat(iii)**3)/(3*10*(1-0.3**2)))**0.25
      
      W=(9.81*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)
      end do

      FSD(n,i,1:nr,iii)=FSD(n-1,i,1:nr,iii) -             &
                           F*FSD(n-1,i,1:nr,iii) + 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(n,i,:,ii)*itd(i,ii)*dx**2/middle_floe_cat**2
      Dave1(ii)=sum(FSD(n,i,:,ii)*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