Blame view

src/initialization.f90 4.99 KB
a081b91c   Jérémy Baudry   new
1
!
c2c1e911   Dany Dumont   Renommage de vari...
2
!____________________________________________________________________
a081b91c   Jérémy Baudry   new
3
!
c2c1e911   Dany Dumont   Renommage de vari...
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
!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

   complex(kind=8), dimension(6) :: poly
   complex(kind=8), dimension(5) :: roots

   logical                       :: param1
   logical                       :: param2

   allocate(Gf(nf))
   allocate(PM(nf))

31a3ff4d   Dany Dumont   Changement de la ...
31
   ! construct time array [in hours]
c2c1e911   Dany Dumont   Renommage de vari...
32
33
34
35
36
   time(1)=0
   do ii=2,nt
      time(ii)=time(ii-1)+dt/3600
   end do

c2c1e911   Dany Dumont   Renommage de vari...
37
  
31a3ff4d   Dany Dumont   Changement de la ...
38
   ! construct space array [in km]
c2c1e911   Dany Dumont   Renommage de vari...
39
40
   x_axis(1)=0
   do ii=2,nx
31a3ff4d   Dany Dumont   Changement de la ...
41
      x_axis(ii)=x_axis(ii-1) + dx/1000
c2c1e911   Dany Dumont   Renommage de vari...
42
   end do
31a3ff4d   Dany Dumont   Changement de la ...
43
44

   write(*,*) '      * dt = ',dt,'s'
37b993a2   Dany Dumont   Nettoyage et reso...
45
   write(*,*) '      * dx = ',dx,'m'
31a3ff4d   Dany Dumont   Changement de la ...
46
47
48
49
50
   if ( disp.eq.1 ) then
      write(*,*) '      * dispersion is on'
   else
      write(*,*) '      * dispersion is off'
   end if
c2c1e911   Dany Dumont   Renommage de vari...
51
 
31a3ff4d   Dany Dumont   Changement de la ...
52
   if ( cont.eq.1 ) then
37b993a2   Dany Dumont   Nettoyage et reso...
53
      write(*,*) '      * wave energy continuously injected'
31a3ff4d   Dany Dumont   Changement de la ...
54
55
56
57
58
   else
      write(*,*) '      * only a single wave energy pulse'
   end if
 

c2c1e911   Dany Dumont   Renommage de vari...
59
60
61
62
63
64
65

!--------------------- INITIAL SPECTRUM -----------------------------

   E(:,:,:)=0d0
   
   if ( init_spec.eq.1 ) then

37b993a2   Dany Dumont   Nettoyage et reso...
66
67
      write(*,*) '      * Jonswap spectrum'

c2c1e911   Dany Dumont   Renommage de vari...
68
69
70
71
72
73
74
75
76
77
      ! JONSWAP spectrum
      do i=1,nf
         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)))
37b993a2   Dany Dumont   Nettoyage et reso...
78
79
      PM=alpha_s*Hs**2*(freq_s**4/freq**5)*                     &
           exp(-beta_s*(freq_s/freq)**4)
c2c1e911   Dany Dumont   Renommage de vari...
80
      Ei=Gf*PM
37b993a2   Dany Dumont   Nettoyage et reso...
81
 
c2c1e911   Dany Dumont   Renommage de vari...
82
   else if ( init_spec.eq.0 ) then
81dede1c   Jérémy Baudry   first commit
83

37b993a2   Dany Dumont   Nettoyage et reso...
84
85
      write(*,*) '      * Bretschneider spectrum'

c2c1e911   Dany Dumont   Renommage de vari...
86
      ! Bretschneider spectrum
37b993a2   Dany Dumont   Nettoyage et reso...
87
      Ei=(1.25*Hs**2*(1/freq)**5)/(8*pi*Tm**4)*                 &
c2c1e911   Dany Dumont   Renommage de vari...
88
           exp(-1.25*((1/freq)/Tm)**4)
a081b91c   Jérémy Baudry   new
89

c2c1e911   Dany Dumont   Renommage de vari...
90
   else
a081b91c   Jérémy Baudry   new
91

37b993a2   Dany Dumont   Nettoyage et reso...
92
93
      Ei=1/(0.01*sqrt(2*pi))*                                   &
           exp(-(omega-2*pi/swell_T)**2/(2*0.01**2))
c2c1e911   Dany Dumont   Renommage de vari...
94
      Ei=(swell_Hs/4d0)**2*Ei/(sum(Ei*domega))
81dede1c   Jérémy Baudry   first commit
95

c2c1e911   Dany Dumont   Renommage de vari...
96
   end if
4508a952   Jérémy Baudry   new version
97

31a3ff4d   Dany Dumont   Changement de la ...
98
   E(1,:,1)=Ei
4508a952   Jérémy Baudry   new version
99

c2c1e911   Dany Dumont   Renommage de vari...
100
!--------------------- ICE TRANSECT ---------------------------------
4508a952   Jérémy Baudry   new version
101

31a3ff4d   Dany Dumont   Changement de la ...
102
   S_ice(1,:,1)=0
c2c1e911   Dany Dumont   Renommage de vari...
103
   x1=floor(X_ice/dx)    !find in which grid bin is the ice edge
37b993a2   Dany Dumont   Nettoyage et reso...
104
   C_ice( 1:x1)=0        !ice concentration is 0 before ice edge
c2c1e911   Dany Dumont   Renommage de vari...
105
   C_ice(x1:nx)=cice     !ice concentration in the transect
4508a952   Jérémy Baudry   new version
106

37b993a2   Dany Dumont   Nettoyage et reso...
107
108
109
   Dave( 1:x1)=0         !initialize mean floe size before ice edge
   Dmax( 1:x1)=0         !initialize max floe size before ice edge
   Dave(x1:nx)=maxfloe   !initialize mean floe size in ice transect
c2c1e911   Dany Dumont   Renommage de vari...
110
   Dmax(x1:nx)=maxfloe   !initialize max floe size in ice transect
4508a952   Jérémy Baudry   new version
111

c2c1e911   Dany Dumont   Renommage de vari...
112
113
   if ( ice_thick.eq.0 ) then
      ! constant ice thickness in the transect
31a3ff4d   Dany Dumont   Changement de la ...
114
      h(x1:nx  )=hice          
c2c1e911   Dany Dumont   Renommage de vari...
115
      h(1 :x1-1)=0d0
4508a952   Jérémy Baudry   new version
116

c2c1e911   Dany Dumont   Renommage de vari...
117
118
119
120
121
122
   else
      ! thickness is a function of distance from ice edge
      h(1:x1)=0d0
      do jj=x1,nx
         h(jj)=hmax*(0.1+0.9*(1-exp(-((real(jj)*dx-X_ice)/Xh))))
      end do
4508a952   Jérémy Baudry   new version
123

c2c1e911   Dany Dumont   Renommage de vari...
124
   end if
4508a952   Jérémy Baudry   new version
125
        
c2c1e911   Dany Dumont   Renommage de vari...
126
!--------------------- FLOE SIZE DISTRIBUTION -----------------------
4e5b4414   Jérémy Baudry   ajout commentaire...
127

37b993a2   Dany Dumont   Nettoyage et reso...
128
   dr=abs(minfloe - maxfloe)/nedge
c2c1e911   Dany Dumont   Renommage de vari...
129
130
   floe_cat(1)=minfloe
   do i=2,nedge
37b993a2   Dany Dumont   Nettoyage et reso...
131
      floe_cat(i)=floe_cat(1)+i*dr
c2c1e911   Dany Dumont   Renommage de vari...
132
   end do
68586e03   Jérémy Baudry   new release
133

c2c1e911   Dany Dumont   Renommage de vari...
134
135
136
   do i=1,nr
      middle_floe_cat(i)=(floe_cat(i)+floe_cat(i+1))*0.5
   end do
68586e03   Jérémy Baudry   new release
137

31a3ff4d   Dany Dumont   Changement de la ...
138
   fsd(:,:,:,:)=0d0
c2c1e911   Dany Dumont   Renommage de vari...
139
   do i=x1,nx
31a3ff4d   Dany Dumont   Changement de la ...
140
      fsd(i,nr,:,1)=C_ice(i)
c2c1e911   Dany Dumont   Renommage de vari...
141
   end do
68586e03   Jérémy Baudry   new release
142

c2c1e911   Dany Dumont   Renommage de vari...
143
   Dave(x1:nx)=middle_floe_cat(nr)
68586e03   Jérémy Baudry   new release
144

c2c1e911   Dany Dumont   Renommage de vari...
145
!--------------------- ICE THICKNESS DISTRIBUTION -------------------
68586e03   Jérémy Baudry   new release
146

37b993a2   Dany Dumont   Nettoyage et reso...
147
   dh=abs(mincat_h - maxcat_h)/nedge_h
c2c1e911   Dany Dumont   Renommage de vari...
148
149
   h_cat(1)=mincat_h
   do i=2,nedge_h
37b993a2   Dany Dumont   Nettoyage et reso...
150
      h_cat(i)=h_cat(1) + i*dh
c2c1e911   Dany Dumont   Renommage de vari...
151
   end do
68586e03   Jérémy Baudry   new release
152

c2c1e911   Dany Dumont   Renommage de vari...
153
   do i=1,nh
37b993a2   Dany Dumont   Nettoyage et reso...
154
      middle_h_cat(i)=0.5*(h_cat(i) + h_cat(i+1))
c2c1e911   Dany Dumont   Renommage de vari...
155
   end do
68586e03   Jérémy Baudry   new release
156

c2c1e911   Dany Dumont   Renommage de vari...
157
   itd(:,:)=0d0
68586e03   Jérémy Baudry   new release
158

c2c1e911   Dany Dumont   Renommage de vari...
159
   if ( itd_scheme.eq.1 ) then
68586e03   Jérémy Baudry   new release
160

c2c1e911   Dany Dumont   Renommage de vari...
161
162
163
164
165
      do i=x1,nx
         itd(i,:)=middle_h_cat/mu_itd**2*                     &
                  exp(-(middle_h_cat-h(i))**2/(2*mu_itd**2))
         itd(i,:)=itd(i,:)/sum(itd(i,:))
      end do
68586e03   Jérémy Baudry   new release
166

c2c1e911   Dany Dumont   Renommage de vari...
167
   else
68586e03   Jérémy Baudry   new release
168

c2c1e911   Dany Dumont   Renommage de vari...
169
170
171
172
173
174
175
      do i=x1,nx
         do j=1,nh
            if (h(i).gt.h_cat(j).and.h(i).le.h_cat(j+1)) then
               itd(i,j)=1
            end if
         end do
      end do
68586e03   Jérémy Baudry   new release
176

c2c1e911   Dany Dumont   Renommage de vari...
177
   end if
68586e03   Jérémy Baudry   new release
178

c2c1e911   Dany Dumont   Renommage de vari...
179
!--------------------- STATISTICS ----------------------------------- 
68586e03   Jérémy Baudry   new release
180

c2c1e911   Dany Dumont   Renommage de vari...
181
   h_sign(1,1)=Hs
68586e03   Jérémy Baudry   new release
182

c2c1e911   Dany Dumont   Renommage de vari...
183
!--------------------- WAVE NUMBER IN ICE ---------------------------
68586e03   Jérémy Baudry   new release
184

c2c1e911   Dany Dumont   Renommage de vari...
185
186
   param1=.false.
   param2=.false.
68586e03   Jérémy Baudry   new release
187
        
c2c1e911   Dany Dumont   Renommage de vari...
188
   open(20,file='input/kice.dat')
68586e03   Jérémy Baudry   new release
189
        
c2c1e911   Dany Dumont   Renommage de vari...
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
   do iii=1,nh
      do ii=1,nf 
         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
208

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