initialization.f90 4.99 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(nf))
   allocate(PM(nf))

   ! construct time array [in hours]
   time(1)=0
   do ii=2,nt
      time(ii)=time(ii-1)+dt/3600
   end do

  
   ! construct space array [in km]
   x_axis(1)=0
   do ii=2,nx
      x_axis(ii)=x_axis(ii-1) + dx/1000
   end do

   write(*,*) '      * dt = ',dt,'s'
   write(*,*) '      * dx = ',dx,'m'
   if ( disp.eq.1 ) then
      write(*,*) '      * dispersion is on'
   else
      write(*,*) '      * dispersion is off'
   end if
 
   if ( cont.eq.1 ) then
      write(*,*) '      * wave energy continuously injected'
   else
      write(*,*) '      * only a single wave energy pulse'
   end if
 


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

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

      write(*,*) '      * Jonswap spectrum'

      ! 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)))
      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

      write(*,*) '      * Bretschneider spectrum'

      ! 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)=Ei

!--------------------- ICE TRANSECT ---------------------------------

   S_ice(1,:,1)=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:nx)=cice     !ice concentration in the transect

   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
   Dmax(x1:nx)=maxfloe   !initialize max floe size in ice transect

   if ( ice_thick.eq.0 ) then
      ! constant ice thickness in the transect
      h(x1:nx  )=hice          
      h(1 :x1-1)=0d0

   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

   end if
        
!--------------------- FLOE SIZE DISTRIBUTION -----------------------

   dr=abs(minfloe - maxfloe)/nedge
   floe_cat(1)=minfloe
   do i=2,nedge
      floe_cat(i)=floe_cat(1)+i*dr
   end do

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

   fsd(:,:,:,:)=0d0
   do i=x1,nx
      fsd(i,nr,:,1)=C_ice(i)
   end do

   Dave(x1:nx)=middle_floe_cat(nr)

!--------------------- ICE THICKNESS DISTRIBUTION -------------------

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

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

   itd(:,:)=0d0

   if ( itd_scheme.eq.1 ) then

      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

   else

      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

   end if

!--------------------- STATISTICS ----------------------------------- 

   h_sign(1,1)=Hs

!--------------------- WAVE NUMBER IN ICE ---------------------------

   param1=.false.
   param2=.false.
        
   open(20,file='input/kice.dat')
        
   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

end subroutine initialization