initialization.f90 4.95 KB
!
!_______________________________________________________________________
!
!               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(nfreq))
                allocate(PM(nfreq))





        !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

        

!_________________________INITIAL SPECTRUM_____________________________

        E(1:nsteps,1:nbin,1:nfreq)=0d0  !INITIALIZE SPECTRUM ARRAY 
        
        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
        
        else if (init_spec.eq.0) then                            !use bretschneider spectrum

        Ei=(1.25*Hs**2*(1/freq)**5)/(8*pi*Tm**4)*exp(-1.25*((1/freq)/Tm)**4)

        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))
        end if



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

!_______________________________________________________________________

!_______________________ICE_TRANSECT____________________________________
S_ice(1,1,1:nfreq)=0
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

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)=maxfloe        !initalize mean floe size in ice transect
Dmax(X1:nbin)=maxfloe        !initialize max floe size in ice transect

if (ice_thick.eq.0) then ! constant ice thickness in the transect

h(X1:nbin)=hice          
h(1:X1-1)=0d0             !0 before ice edge

else                    ! thickness is a function of distance from ice edge

h(1:X1)=0d0
        do jj=X1,nbin

                h(jj)=hmax*(0.1+0.9*(1-exp(-((real(jj)*dx-X_ice)/Xh))))

        end do

end if
        
!_____________________FSD_____________________________________________
        res=abs(minfloe-maxfloe)/nedge
floe_cat(1)=minfloe
do i=2,nedge
floe_cat(i)=floe_cat(1)+i*res
end do

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

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

        








end subroutine initialization