ice_fracture.f90 5.49 KB
Newer Older
1 2 3 4 5 6 7 8 9
!--------------------------------------------------------------------
!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.
Jérémy Baudry's avatar
Jérémy Baudry committed
10

11
subroutine ice_fracture
Jérémy Baudry's avatar
Jérémy Baudry committed
12

13
   use parameters
Jérémy Baudry's avatar
Jérémy Baudry committed
14

15 16
   !local parameters
   implicit none
Jérémy Baudry's avatar
Jérémy Baudry committed
17
        
18
   !integer                         :: Nband=500
19 20 21 22
   integer                         :: freq1
   integer                         :: freq2
   integer                         :: ninterp=1

23
   !double precision                :: TT=13
24 25 26 27
   double precision                :: m_0
   double precision                :: m_0_stress
   double precision                :: m_0_amp
   double precision                :: m_2_amp
28 29
   !double precision                :: Tw
   !double precision                :: kw
30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
   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
80 81 82 83 84 85 86
      !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
87
      wavelength=2*pi/kice(iii,:)
Jérémy Baudry's avatar
Jérémy Baudry committed
88

89
      strain=middle_h_cat(iii)*0.5*kice(iii,:)**2*W
Jérémy Baudry's avatar
Jérémy Baudry committed
90

91 92
      freq1=nf
      freq2=nf
Jérémy Baudry's avatar
Jérémy Baudry committed
93

94
      do ii=1,nr
Jérémy Baudry's avatar
Jérémy Baudry committed
95
        
96 97
         if ( 2d0*floe_cat(ii+1).lt.wavelength(nf) .or.         &
              2d0*floe_cat(ii+1).gt.wavelength( 1) ) then
98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123
            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      
Jérémy Baudry's avatar
Jérémy Baudry committed
124
     
125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142
      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
143 144
         !Q2(jj)=sum(Q(jj,:)*fsd(n-1,i,:,iii)*F)
         Q2(jj)=sum(Q(jj,:)*fsd(i,:,iii,n-1)*F)
145 146
      end do

147 148 149 150
      !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
Jérémy Baudry's avatar
Jérémy Baudry committed
151
        
152 153 154 155
      !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')
Jérémy Baudry's avatar
Jérémy Baudry committed
156
        
157 158 159 160 161 162 163 164 165 166
      !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
167 168
      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)
169 170 171 172 173 174 175 176 177 178 179 180
   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