Blame view

src/initialization.f90 3.05 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
20
21
22
!
!_______________________________________________________________________
!
!               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

                allocate(Gf(nfreq))
                allocate(PM(nfreq))
81dede1c   Jérémy Baudry   first commit
23

81dede1c   Jérémy Baudry   first commit
24

81dede1c   Jérémy Baudry   first commit
25

4e5b4414   Jérémy Baudry   ajout commentaire...
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41


        !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
42
!_________________________INITIAL SPECTRUM_____________________________
a081b91c   Jérémy Baudry   new
43
44

        E(1:nsteps,1:nbin,1:nfreq)=0d0  !INITIALIZE SPECTRUM ARRAY 
4508a952   Jérémy Baudry   new version
45
        
a081b91c   Jérémy Baudry   new
46
47
48
49
50
51
52
53
54
55
56
57
58
59
        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
60
        
a081b91c   Jérémy Baudry   new
61
        else                            !use bretschneider spectrum
4508a952   Jérémy Baudry   new version
62
63
64

        Ei=(1.25*Hs**2*(1/freq)**5)/(8*pi*Tm**4)*exp(-1.25*((1/freq)/Tm)**4)
        end if
81dede1c   Jérémy Baudry   first commit
65

81dede1c   Jérémy Baudry   first commit
66

a081b91c   Jérémy Baudry   new
67
68
69

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

81dede1c   Jérémy Baudry   first commit
70
71
!_______________________________________________________________________

4508a952   Jérémy Baudry   new version
72
!_______________________ICE_TRANSECT____________________________________
4e5b4414   Jérémy Baudry   ajout commentaire...
73
S_ice(1,1,1:nfreq)=0
a081b91c   Jérémy Baudry   new
74
75
76
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
77

a081b91c   Jérémy Baudry   new
78
79
80
81
Dave(1:X1)=0            !initalize mean floe size before ice edge
Dmax(1:X1)=0            !initialize max floe size before ice edge
Dave(X1:nbin)=D0        !initalize mean floe size in ice transect
Dmax(X1:nbin)=D0        !initialize max floe size in ice transect
4508a952   Jérémy Baudry   new version
82

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

a081b91c   Jérémy Baudry   new
85
86
h(X1:nbin)=hice          
h(1:X1)=0d0             !0 before ice edge
4508a952   Jérémy Baudry   new version
87

a081b91c   Jérémy Baudry   new
88
else                    ! thickness is a function of distance from ice edge
4508a952   Jérémy Baudry   new version
89
90

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

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

a081b91c   Jérémy Baudry   new
95
        end do
4508a952   Jérémy Baudry   new version
96
97
98
99

end if
        
!_____________________FSD_____________________________________________
4e5b4414   Jérémy Baudry   ajout commentaire...
100
101
102
103
104
105
106
107
108
109
110
        res=abs(minfloe-maxfloe)/nedge

floe_cat(1)=maxfloe
do i=2,nedge
floe_cat(i)=floe_cat(1)-i*res
end do

do i=1,nedge-1
middle_floe_cat(i)=(floe_cat(i)+floe_cat(i+1))*0.5
end do

4508a952   Jérémy Baudry   new version
111
112
113

FSD(1,1:nbin,1)=1d0

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