Blame view

src/initialization.f90 4.95 KB
a081b91c   Jérémy Baudry   new
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
!
!_______________________________________________________________________
!
!               DESCRIPTION:
!               This is the initialization routine which  contains the 
!               initial spectrum construction, set initial values for arrays
!               and construct the ice transect.
!_______________________________________________________________________

                !INTERFACE:
                subroutine initialization

                !MODULE USES:
                use parameters

                !LOCAL PARAMETERS AND VARIABLES
                implicit none
                double precision, allocatable ::Gf(:),PM(:)
                integer ::X1
68586e03   Jérémy Baudry   new release
20
21
22
23
                complex(kind=8), dimension(6)  :: poly
                complex(kind=8), dimension(5)  :: roots
                logical                        :: param1
                logical                        :: param2
a081b91c   Jérémy Baudry   new
24
25
                allocate(Gf(nfreq))
                allocate(PM(nfreq))
81dede1c   Jérémy Baudry   first commit
26

81dede1c   Jérémy Baudry   first commit
27

81dede1c   Jérémy Baudry   first commit
28

4e5b4414   Jérémy Baudry   ajout commentaire...
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44


        !construct time array
         time(1)=0
        do ii=2,nsteps
        time(ii)=time(ii-1)+dt/60
        end do
        
        !construct space array
        x_axis(1)=0
        do ii=2,nbin
        x_axis(ii)=x_axis(ii-1)+dx/1000
        end do

        

81dede1c   Jérémy Baudry   first commit
45
!_________________________INITIAL SPECTRUM_____________________________
a081b91c   Jérémy Baudry   new
46
47

        E(1:nsteps,1:nbin,1:nfreq)=0d0  !INITIALIZE SPECTRUM ARRAY 
4508a952   Jérémy Baudry   new version
48
        
a081b91c   Jérémy Baudry   new
49
50
51
52
53
54
55
56
57
58
59
60
61
62
        if(init_spec.eq.1) then         !use JONSWAP spectrum

        do i=1,nfreq 
                if (freq(i).le.freq_s) then
                        sigma_s(i)=0.07
                else
                        sigma_s(i)=0.09
                end if
        end do

        Gf=gamma_s**(exp((-(freq-freq_s)**2)/(2*sigma_s**2*freq_s**2)))
        PM=alpha_s*Hs**2*(freq_s**4/freq**5)*exp(-beta_s*(freq_s/freq)**4)

        Ei=Gf*PM
4508a952   Jérémy Baudry   new version
63
        
68586e03   Jérémy Baudry   new release
64
        else if (init_spec.eq.0) then                            !use bretschneider spectrum
4508a952   Jérémy Baudry   new version
65
66

        Ei=(1.25*Hs**2*(1/freq)**5)/(8*pi*Tm**4)*exp(-1.25*((1/freq)/Tm)**4)
68586e03   Jérémy Baudry   new release
67
68
69
70

        else
        Ei=1/(0.01*sqrt(2*pi))*exp(-(omega-2*pi/swell_T)**2/(2*0.01**2))
        Ei=(swell_Hs/4d0)**2*Ei/(sum(Ei*domega))
4508a952   Jérémy Baudry   new version
71
        end if
81dede1c   Jérémy Baudry   first commit
72

81dede1c   Jérémy Baudry   first commit
73

a081b91c   Jérémy Baudry   new
74
75
76

        E(1,1,1:nfreq)=Ei               !initial spectrum

81dede1c   Jérémy Baudry   first commit
77
78
!_______________________________________________________________________

4508a952   Jérémy Baudry   new version
79
!_______________________ICE_TRANSECT____________________________________
4e5b4414   Jérémy Baudry   ajout commentaire...
80
S_ice(1,1,1:nfreq)=0
a081b91c   Jérémy Baudry   new
81
82
83
X1=floor(X_ice/dx)      !find in which grid bin is the ice edge
C_ice(1:X1)=0           !ice concentration is 0 before ice edge!
C_ice(X1:nbin)=cice     !ice concentration in the transect
4508a952   Jérémy Baudry   new version
84

a081b91c   Jérémy Baudry   new
85
86
Dave(1:X1)=0            !initalize mean floe size before ice edge
Dmax(1:X1)=0            !initialize max floe size before ice edge
68586e03   Jérémy Baudry   new release
87
88
Dave(X1:nbin)=maxfloe        !initalize mean floe size in ice transect
Dmax(X1:nbin)=maxfloe        !initialize max floe size in ice transect
4508a952   Jérémy Baudry   new version
89

a081b91c   Jérémy Baudry   new
90
if (ice_thick.eq.0) then ! constant ice thickness in the transect
4508a952   Jérémy Baudry   new version
91

a081b91c   Jérémy Baudry   new
92
h(X1:nbin)=hice          
68586e03   Jérémy Baudry   new release
93
h(1:X1-1)=0d0             !0 before ice edge
4508a952   Jérémy Baudry   new version
94

a081b91c   Jérémy Baudry   new
95
else                    ! thickness is a function of distance from ice edge
4508a952   Jérémy Baudry   new version
96
97

h(1:X1)=0d0
a081b91c   Jérémy Baudry   new
98
        do jj=X1,nbin
4508a952   Jérémy Baudry   new version
99

a081b91c   Jérémy Baudry   new
100
                h(jj)=hmax*(0.1+0.9*(1-exp(-((real(jj)*dx-X_ice)/Xh))))
4508a952   Jérémy Baudry   new version
101

a081b91c   Jérémy Baudry   new
102
        end do
4508a952   Jérémy Baudry   new version
103
104
105
106

end if
        
!_____________________FSD_____________________________________________
4e5b4414   Jérémy Baudry   ajout commentaire...
107
        res=abs(minfloe-maxfloe)/nedge
68586e03   Jérémy Baudry   new release
108
floe_cat(1)=minfloe
4e5b4414   Jérémy Baudry   ajout commentaire...
109
do i=2,nedge
68586e03   Jérémy Baudry   new release
110
floe_cat(i)=floe_cat(1)+i*res
4e5b4414   Jérémy Baudry   ajout commentaire...
111
112
end do

68586e03   Jérémy Baudry   new release
113
do i=1,nbcat
4e5b4414   Jérémy Baudry   ajout commentaire...
114
115
116
middle_floe_cat(i)=(floe_cat(i)+floe_cat(i+1))*0.5
end do

68586e03   Jérémy Baudry   new release
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
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
FSD(:,:,:,:)=0d0
do i=X1,Nbin
FSD(1,i,nbcat,:)=C_ice(i)
end do

Dave(X1:nbin)=middle_floe_cat(nbcat)        !initialize max floe size in ice transect



!_____________________IDT_____________________________________________
        resh=abs(mincat_h-maxcat_h)/(nedge_h)
h_cat(1)=mincat_h
do i=2,nedge_h
h_cat(i)=h_cat(1)+i*resh
end do

do i=1,nbcat_h
middle_h_cat(i)=(h_cat(i)+h_cat(i+1))*0.5
end do

IDT(:,:)=0d0

if(IDT_scheme.eq.1) then

do i=X1,nbin
IDT(i,:)=middle_h_cat/mu_IDT**2*exp(-(middle_h_cat-h(i))**2/(2*mu_IDT**2))
IDT(i,:)=IDT(i,:)/sum(IDT(i,:))
end do

else

        do i=X1,nbin
        do j=1,nbcat_h
        if (h(i).gt.h_cat(j).and.h(i).le.h_cat(j+1)) then
        IDT(i,j)=1
        end if
        end do
        end do

end if





!_________________STATISTICS_____________________________________________

h_sign(1,1)=Hs
!_________________wave number in ice ____________________________________

        param1=.false.
        param2=.false.
        
           open(20,file='output/kice.dat')
        
        do iii=1,nbcat_h
        do  ii=1,nfreq 
        poly(1)=-922*omega(ii)**2
        poly(2)=922*9.81-0.9*middle_h_cat(iii)*omega(ii)**2
        poly(3)=0
        poly(4)=0
        poly(5)=0
        poly(6)=(5.5e9*middle_h_cat(iii)**3)/(12*(1-0.3**2))

        call cmplx_roots_gen(roots,poly,5,param1,param2)
        do jj=1,5
        if (real(roots(jj)).gt.0 .and. aimag(roots(jj)).eq.0) then
        kice(iii,ii)=real(roots(jj))
        end if
        end do
        end do
        write(20,*)kice(iii,:)
        end do

        






4508a952   Jérémy Baudry   new version
198

4508a952   Jérémy Baudry   new version
199

81dede1c   Jérémy Baudry   first commit
200
end subroutine initialization