Blame view

src/ice_fracture.f90 5.11 KB
68586e03   Jérémy Baudry   new release
1
2
3
4
5
6
7
8
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
74
75
76
77
78
79
80
81
82
83
84
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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
      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