initialization.f90 4.95 KB
Newer Older
Jérémy Baudry's avatar
new  
Jérémy Baudry committed
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
Jérémy Baudry's avatar
Jérémy Baudry committed
20 21 22 23
                complex(kind=8), dimension(6)  :: poly
                complex(kind=8), dimension(5)  :: roots
                logical                        :: param1
                logical                        :: param2
Jérémy Baudry's avatar
new  
Jérémy Baudry committed
24 25
                allocate(Gf(nfreq))
                allocate(PM(nfreq))
Jérémy Baudry's avatar
Jérémy Baudry committed
26 27 28



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

        

Jérémy Baudry's avatar
Jérémy Baudry committed
45
!_________________________INITIAL SPECTRUM_____________________________
Jérémy Baudry's avatar
new  
Jérémy Baudry committed
46 47

        E(1:nsteps,1:nbin,1:nfreq)=0d0  !INITIALIZE SPECTRUM ARRAY 
Jérémy Baudry's avatar
Jérémy Baudry committed
48
        
Jérémy Baudry's avatar
new  
Jérémy Baudry committed
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
Jérémy Baudry's avatar
Jérémy Baudry committed
63
        
Jérémy Baudry's avatar
Jérémy Baudry committed
64
        else if (init_spec.eq.0) then                            !use bretschneider spectrum
Jérémy Baudry's avatar
Jérémy Baudry committed
65 66

        Ei=(1.25*Hs**2*(1/freq)**5)/(8*pi*Tm**4)*exp(-1.25*((1/freq)/Tm)**4)
Jérémy Baudry's avatar
Jérémy Baudry committed
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))
Jérémy Baudry's avatar
Jérémy Baudry committed
71
        end if
Jérémy Baudry's avatar
Jérémy Baudry committed
72 73


Jérémy Baudry's avatar
new  
Jérémy Baudry committed
74 75 76

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

Jérémy Baudry's avatar
Jérémy Baudry committed
77 78
!_______________________________________________________________________

Jérémy Baudry's avatar
Jérémy Baudry committed
79
!_______________________ICE_TRANSECT____________________________________
80
S_ice(1,1,1:nfreq)=0
Jérémy Baudry's avatar
new  
Jérémy Baudry committed
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
Jérémy Baudry's avatar
Jérémy Baudry committed
84

Jérémy Baudry's avatar
new  
Jérémy Baudry committed
85 86
Dave(1:X1)=0            !initalize mean floe size before ice edge
Dmax(1:X1)=0            !initialize max floe size before ice edge
Jérémy Baudry's avatar
Jérémy Baudry committed
87 88
Dave(X1:nbin)=maxfloe        !initalize mean floe size in ice transect
Dmax(X1:nbin)=maxfloe        !initialize max floe size in ice transect
Jérémy Baudry's avatar
Jérémy Baudry committed
89

Jérémy Baudry's avatar
new  
Jérémy Baudry committed
90
if (ice_thick.eq.0) then ! constant ice thickness in the transect
Jérémy Baudry's avatar
Jérémy Baudry committed
91

Jérémy Baudry's avatar
new  
Jérémy Baudry committed
92
h(X1:nbin)=hice          
Jérémy Baudry's avatar
Jérémy Baudry committed
93
h(1:X1-1)=0d0             !0 before ice edge
Jérémy Baudry's avatar
Jérémy Baudry committed
94

Jérémy Baudry's avatar
new  
Jérémy Baudry committed
95
else                    ! thickness is a function of distance from ice edge
Jérémy Baudry's avatar
Jérémy Baudry committed
96 97

h(1:X1)=0d0
Jérémy Baudry's avatar
new  
Jérémy Baudry committed
98
        do jj=X1,nbin
Jérémy Baudry's avatar
Jérémy Baudry committed
99

Jérémy Baudry's avatar
new  
Jérémy Baudry committed
100
                h(jj)=hmax*(0.1+0.9*(1-exp(-((real(jj)*dx-X_ice)/Xh))))
Jérémy Baudry's avatar
Jérémy Baudry committed
101

Jérémy Baudry's avatar
new  
Jérémy Baudry committed
102
        end do
Jérémy Baudry's avatar
Jérémy Baudry committed
103 104 105 106

end if
        
!_____________________FSD_____________________________________________
107
        res=abs(minfloe-maxfloe)/nedge
Jérémy Baudry's avatar
Jérémy Baudry committed
108
floe_cat(1)=minfloe
109
do i=2,nedge
Jérémy Baudry's avatar
Jérémy Baudry committed
110
floe_cat(i)=floe_cat(1)+i*res
111 112
end do

Jérémy Baudry's avatar
Jérémy Baudry committed
113
do i=1,nbcat
114 115 116
middle_floe_cat(i)=(floe_cat(i)+floe_cat(i+1))*0.5
end do

Jérémy Baudry's avatar
Jérémy Baudry committed
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

        






Jérémy Baudry's avatar
Jérémy Baudry committed
198 199


Jérémy Baudry's avatar
Jérémy Baudry committed
200
end subroutine initialization