ice_fracture.f90 5.11 KB
      subroutine ice_fracture

        use parameters

        !local parameters
        implicit none
        
        double precision                :: Y=5.5e9
        double precision                :: TT=13
        double precision                :: poisson=0.3
        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(:)
        integer                         :: Nband=500
        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,allocatable    :: wavelength(:)
        double precision,allocatable    :: S_strain(:)
        double precision, allocatable   :: S_stress(:)
        integer                         :: freq1
        integer                         :: freq2
        double precision, allocatable   :: k(:)
        double precision                :: Lmin
        double precision, allocatable   :: Dave1(:)
        integer                         ::ninterp=1
        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(P_cat(nbcat))
        allocate(F(nbcat))
        allocate(Q(nbcat,nbcat))
        allocate(E_A(nfreq))
        allocate(Q2(nbcat))
        allocate(strain(nfreq))
        allocate(stress(nfreq))
        allocate(wavelength(nfreq))
        allocate(S_strain(nbcat))
        allocate(k(nfreq))
        allocate(S_stress(Nband))
        allocate(Dave1(nbcat_h))
        allocate(xinterp(ninterp))
        allocate(Pinterp(ninterp))
        allocate(nfloe(nbcat,nbcat_h))
        allocate(nfloe_tot(nbcat))
        allocate(P2(Nband))
        allocate(W(nfreq))
       
      

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


        do iii=1,nbcat_h

          Lmin=pi*0.5*((5e6*middle_h_cat(iii)**3)/(3*10*(1-0.3**2)))**0.25!Lmin according Mellor 1986

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

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

        freq1=nfreq
        freq2=nfreq

        do ii=1,nbcat
        
        if(2d0*floe_cat(ii+1).lt.wavelength(nfreq).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,nbcat
        F(ii)=sum(middle_floe_cat(1:ii)*P_cat(1:ii)/dx)
        end do

      

        Q(:,:)=0
        !this is the redistribution function
        Q(:,1)=0 ! the smallest size category cannot redistribute
        do ii=2,nbcat
        !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,nbcat

        Q2(jj)=sum(Q(jj,:)*FSD(n-1,i,:,iii)*F)
        end do

        FSD(n,i,1:nbcat,iii)=FSD(n-1,i,1:nbcat,iii)-F*FSD(n-1,i,1:nbcat,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,nbcat
        !write(32,*)Q(:,ii)
        !end do
        !end if
        end do
        !compute the new average floe diameter <D>
        do ii=1,nbcat_h
        nfloe(:,ii)=FSD(n,i,:,ii)*IDT(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,nbcat
        nfloe_tot(ii)=sum(nfloe(ii,:))
        end do
        Dave(i)=sum(middle_floe_cat*nfloe_tot)
        !Dave(i)=sum(IDT(i,:)*Dave1)
       end subroutine