Blame view

src/ice_fracture.f90 4.84 KB
68586e03   Jérémy Baudry   new release
1

c2c1e911   Dany Dumont   Renommage de vari...
2
subroutine ice_fracture
68586e03   Jérémy Baudry   new release
3

c2c1e911   Dany Dumont   Renommage de vari...
4
   use parameters
68586e03   Jérémy Baudry   new release
5

c2c1e911   Dany Dumont   Renommage de vari...
6
7
   !local parameters
   implicit none
68586e03   Jérémy Baudry   new release
8
        
c2c1e911   Dany Dumont   Renommage de vari...
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
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
   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
68586e03   Jérémy Baudry   new release
74
      
c2c1e911   Dany Dumont   Renommage de vari...
75
76
      W=(9.81*kice(iii,:))/omega**2
      wavelength=2*pi/kice(iii,:)
68586e03   Jérémy Baudry   new release
77

c2c1e911   Dany Dumont   Renommage de vari...
78
      strain=middle_h_cat(iii)*0.5*kice(iii,:)**2*W
68586e03   Jérémy Baudry   new release
79

c2c1e911   Dany Dumont   Renommage de vari...
80
81
      freq1=nf
      freq2=nf
68586e03   Jérémy Baudry   new release
82

c2c1e911   Dany Dumont   Renommage de vari...
83
      do ii=1,nr
68586e03   Jérémy Baudry   new release
84
        
c2c1e911   Dany Dumont   Renommage de vari...
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
         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      
68586e03   Jérémy Baudry   new release
113
     
c2c1e911   Dany Dumont   Renommage de vari...
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
      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
68586e03   Jérémy Baudry   new release
137
        
c2c1e911   Dany Dumont   Renommage de vari...
138
139
140
141
      !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')
68586e03   Jérémy Baudry   new release
142
        
c2c1e911   Dany Dumont   Renommage de vari...
143
144
145
146
147
148
149
150
151
152
153
154
155
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
      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