first commit

nml/parameter.nml

+
+&waves_parameters
+
+Tm	=6
+Hs	=1
+disp	=0
+
+
+
+/
+
+&model_parameter
+
+nbin	=10
+dx	=5000
+Cfl	=0.7
+/
+
+
+&spectrum_parameters
+
+
+nfreq	=60
+Tmin	=2.5
+Tmax	=20
+alpha_s =0.2044
+beta_s	=1.2500
+gamma_s	=3.3
+
+
+/
+
+&ice_parameters
+cice	=0.5
+hice	=1
+D0	=500
+gam	=1.5
+Dmin	=20
+
+/

nml/parameter.nml~

+
+&waves_parameters
+
+Tm	=6
+Hs	=1
+disp	=0
+
+
+
+/
+
+&model_parameter
+
+nbin	=10
+dx	=5000
+Cfl	=0.1
+/
+
+
+&spectrum_parameters
+
+
+nfreq	=60
+Tmin	=2.5
+Tmax	=20
+alpha_s =0.2044
+beta_s	=1.2500
+gamma_s	=3.3
+
+
+/
+
+&ice_parameters
+cice	=0.5
+hice	=1
+D0	=500
+gam	=1.5
+Dmin	=20
+
+/

output/data_treatment.m

+% clear all;
+close all;
+
+E=load('output/Energy_spectrum.dat');
+subplot(1,2,1)
+cmap=jet(size(E,1));
+for i=1:size(E,1)
+    plot(E(i,:),'color',cmap(i,:))
+    hold on
+    
+  
+end
+
+FSD=load('output/floe_size.dat');
+subplot(1,2,2)
+plot(FSD(:,1))
+hold on
+plot(FSD(:,2),'r')

src/advection.f90

+subroutine advection
+
+use parameters
+
+implicit none
+
+	double precision, allocatable :: beta_1(:),beta_2(:),beta_3(:)
+
+
+	allocate(beta_1(nfreq))
+	allocate(beta_2(nfreq))
+	allocate(beta_3(nfreq))
+
+	beta_1=CN/2*(CN+1)
+	beta_2=1-CN**2
+	beta_3=CN/2*(CN-1)
+
+	
+	E(n,1,1:nfreq)=beta_2*E(n-1,1,1:nfreq)+beta_3*E(n-1,2,1:nfreq)
+
+	if (i.le.nbin) then
+	E(n,i,1:nfreq)=beta_1*E(n-1,i-1,1:nfreq)+beta_2*E(n-1,i,1:nfreq)+beta_3*E(n-1,i+1,1:nfreq)
+	else
+		
+	E(n,i,1:nfreq)=beta_1*E(n-1,i-1,1:nfreq)!+beta_2*E(n-1,i,1:nfreq)
+	end if
+	
+	
+	do j=1,nbin
+	do ii=1,nfreq
+	if (E(n,j,ii).lt.0) then
+		E(n,j,ii)=0
+	end if
+	end do
+	end do
+	
+
+end subroutine advection

src/attenuation.f90

+subroutine attenuation
+use parameters
+
+implicit none
+
+	double precision :: q11,q12,q13,q14,q15,q21,q22,q23,q24,q25
+	double precision :: q31,q32,q33,q34,q35
+	double precision, allocatable ::p1(:),p2(:),p3(:),att(:)
+
+	allocate(p1(nfreq))
+	allocate(p2(nfreq))
+	allocate(p3(nfreq))
+	allocate(att(nfreq))
+
+
+q11 = -0.00000777
+q12 =  0.00032080
+q13 = -0.00437542
+q14 =  0.02047559
+q15 = -0.01356537
+
+q21 =  0.00003635
+q22 = -0.00153484
+q23 =  0.02121709
+q24 = -0.09289399
+q25 = -0.03693082
+
+q31 = -0.00004509
+q32 =  0.00214484
+q33 = -0.03663425
+q34 =  0.26065369
+q35 = -0.62474085
+
+p1 = q11*T**4 + q12*T**3 + q13*T**2 + q14*T + q15
+p2 = q21*T**4 + q22*T**3 + q23*T**2 + q24*T + q25
+p3 = q31*T**4 + q32*T**3 + q33*T**2 + q34*T + q35
+
+alpha(i,1:nfreq)=-1*(p1*h(i)**2 + p2*h(i) + p3)
+
+do j=1,nfreq
+if (alpha(i,j).lt.0)then
+alpha(i,j)=0
+end if
+end do
+
+att=C_ice(i)*alpha(i,1:nfreq)/Dave(i)
+
+
+S_ice(i,1:nfreq)=-att*E(n,i,1:nfreq)
+
+E(n,i,1:nfreq)=E(n,i,1:nfreq)*exp(-att*Cg*dt)
+
+
+
+
+end subroutine attenuation

src/floe_breaking.f90

+subroutine floe_breaking
+
+use parameters
+
+implicit none
+	double precision, allocatable :: WN_ice(:),fnt(:),k(:)
+	double precision ::rho=1025,d,F,Y=5.5,v=0.3
+	integer :: acc=1000
+	double precision, allocatable :: W(:)
+	double precision, allocatable :: fnt_d(:)
+	
+	double precision, allocatable :: S(:)
+	double precision :: SS
+	double precision :: WN_ice_d
+	double precision :: RS=7.05e-05
+	double precision :: m0_s
+	double precision :: m0_a
+	double precision :: m2_a
+	double precision :: Tw
+	double precision :: om_d
+	double precision :: wl_w
+
+	allocate(fnt(acc+1))
+	allocate(k(acc+1))
+	allocate(WN_ice(nfreq))
+	allocate(W(nfreq))
+	allocate(fnt_d(acc+1))
+	allocate(S(nfreq))
+
+
+	F=Y*h(i)**3/12*(1-v**2)
+	d=0.9*h(i)
+
+	jj=1
+	do ii=0,acc
+	k(jj)=real(ii)*2/real(acc)
+	jj=jj+1
+	end do
+	
+	
+	do j=1,nfreq
+	fnt=F*k**5+rho*(g-d*omega(j)**2)*k-rho*omega(j)**2
+	WN_ice(j)=k(minloc(abs(fnt),DIM=1))
+	
+	end do
+
+ 
+	W=g*WN_ice/omega**2
+	S=(h(i)/2)*WN_ice**2*W
+	m0_s=sum(E(n,i,1:nfreq)*S**2)*domega
+	SS=2*sqrt(m0_s)
+
+	if (SS.gt.RS) then
+
+	m0_a=sum(E(n,i,1:nfreq)*W**2)*domega
+	m2_a=sum(omega**2*E(n,i,1:nfreq)*W**2)*domega
+	Tw=2*pi*sqrt(m0_a/m2_a)
+	om_d=2*pi/Tw
+
+	fnt_d=F*k**5+rho*(g-d*om_d**2)*k-rho*om_d**2
+	WN_ice_d=k(minloc(abs(fnt_d),DIM=1))
+
+	wl_w=2*pi/WN_ice_d
+
+	Dmax(i)=max(Dmin,min(wl_w/2,Dmax(i)))
+	
+
+	end if
+
+
+
+
+end subroutine floe_breaking

src/fsd_build.f90

+subroutine fsd_build
+use parameters
+
+	implicit none
+	double precision :: coeff
+
+	if (Dmax(i).eq.Dmin) then
+	Dave(i)=Dmin
+	elseif (Dmax(i).eq.D0) then
+	Dave(i)=Dmax(i)
+	else
+	coeff=1/((1/(1-gam))*(Dmax(i)**(1-gam)-Dmin**(1-gam)))
+	Dave(i)=coeff*(1/(2-gam))*(Dmax(i)**(2-gam)-Dmin**(2-gam))
+	end if
+
+
+end subroutine fsd_build

src/initialization.f90

+subroutine initialization
+use parameters
+
+!local parameters
+implicit none
+	double precision, allocatable ::Gf(:),PM(:)
+
+	allocate(Gf(nfreq))
+	allocate(PM(nfreq))
+	
+	
+
+!_________________________INITIAL SPECTRUM_____________________________
+	E(1,1:nbin,1:nfreq)=0
+
+	!build 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
+
+	
+
+	E(1,1,1:nfreq)=Ei
+	
+!_______________________________________________________________________
+
+end subroutine initialization

src/main.f90

+PROGRAM WIM2
+
+
+	use parameters
+	call read_namelist
+	call array_allocation
+	call initialization
+	
+	
+
+
+!______________________________________________________
+
+
+
+
+
+!_______________________TIME LOOP______________________
+	do n=2,nsteps
+	
+	call progress(n,nsteps)
+	
+		do i=2,nbin
+		call advection		
+		!call attenuation
+		!call floe_breaking
+		!call fsd_build
+		end do
+	end do
+	!______________________OUTPUTS_________________________
+	root='output/'
+	spectrum=trim(root)//'Energy_spectrum.dat'
+	floe_size=trim(root)//'floe_size.dat'
+
+
+	open(10,file=spectrum)
+	open(11,file=floe_size)
+
+	!call write_output
+	do i=1,nbin
+	!do j=1,nsteps
+	write(10,*)E(nsteps,i,1:nfreq)
+	!end do
+	write(11,*)Dave(i),Dmax(i)
+	end do
+
+	close(10)
+	close(11)
+
+
+
+contains
+
+subroutine progress(j,jmax)
+  implicit none
+  integer::j,k,jmax
+  character(len=15)::bar="processing ???%"
+  write(unit=bar(12:14),fmt="(i3)")ceiling((real(j)/real(jmax))*100)
+  write(unit=6,fmt="(a1,a17)",advance="no") char(13), bar
+  if (ceiling((real(j)/real(jmax))*100).eq.100) then
+  write(unit=6,fmt=*)
+  write(*,*)'simulation completed!'
+  
+  endif
+end subroutine progress
+
+END PROGRAM WIM2

src/makefile

+COMPILER= gfortran
+
+OPTION= -O3 
+
+
+NETCDFinc= -I/usr/local/include 
+NETCDFLIB= -L/usr/local/lib -lnetcdf
+
+OBJ= *.o
+
+SRC= *.f90
+
+MOD= parameters.mod 
+
+EXEC= WIM2
+
+
+
+WIM2: $(OBJ)
+	$(COMPILER) $(OPTION) $(NETCDFinc) -o $(EXEC) $(OBJ) $(NETCDFLIB) 
+	mv WIM2 ../
+	
+
+$(MOD): parameters.f90 
+	$(COMPILER)  -c parameters.f90
+
+
+$(OBJ): $(SRC) $(MOD)
+	$(COMPILER) $(NETCDFinc) -c $(SRC) $(NETCDFLIB) 
+
+
+clean:
+	rm -f *.o *.mod 
+
+mrproper:
+
+	rm -f *.o *.mod ../WIM2

src/makefile~

+COMPILER= ifort
+
+OPTION= -O3 
+
+
+NETCDFinc= -I/usr/local/include 
+NETCDFLIB= -L/usr/local/lib -lnetcdf
+
+OBJ= *.o
+
+SRC= *.f90
+
+MOD= parameters.mod 
+
+EXEC= WIM2
+
+
+
+WIM2: $(OBJ)
+	$(COMPILER) $(OPTION) $(NETCDFinc) -o $(EXEC) $(OBJ) $(NETCDFLIB) 
+	mv WIM2 ../
+	
+
+$(MOD): parameters.f90 
+	$(COMPILER)  -c parameters.f90
+
+
+$(OBJ): $(SRC) $(MOD)
+	$(COMPILER) $(NETCDFinc) -c $(SRC) $(NETCDFLIB) 
+
+
+clean:
+	rm -f *.o *.mod 
+
+mrproper:
+
+	rm -f *.o *.mod ../WIM2

src/parameters.f90

+module parameters
+
+implicit none
+
+!__________________name of files__________________
+	character(len=100) ::root,spectrum,floe_size
+
+!______________dummys_____________________________
+	integer	:: i,ii,n,j,jj
+!______________global parameters__________________
+
+	double precision	:: g=9.81,pi=3.1416
+!____________Spectrum_____________________________
+	double precision, allocatable	:: E(:,:,:),Ei(:),T(:)
+	double precision, allocatable	:: omega(:),freq(:),sigma_s(:)
+	integer	:: nbin=10,nfreq=60
+	double precision :: alpha_s=8.1e-3,beta_s=1.25
+	double precision :: Tmin=2.5,Tmax=20,Tm=6,Hs=1
+	double precision :: gamma_s=3.3,freq_s,domega
+!____________waves_________________________________
+	double precision, allocatable	::wl(:),Cp(:),Cg(:),CN(:)
+	integer				::disp=0
+!____________grid parameters________________________
+	double precision :: dt,dx=500
+	integer		 :: nsteps
+	double precision :: Cfl=1
+!_____________ICE_____________________________________
+	double precision, allocatable :: alpha(:,:)
+	double precision, allocatable :: h(:)
+	double precision, allocatable :: C_ice(:)
+	double precision, allocatable :: S_ice(:,:)
+	double precision :: cice=1
+	double precision :: hice=1
+	double precision, allocatable :: Dave(:)
+	double precision, allocatable :: Dmax(:)
+	double precision :: D0=200
+	double precision :: gam=2
+	double precision :: Dmin=20
+	
+
+	
+
+
+contains
+
+
+subroutine read_namelist
+
+	namelist /spectrum_parameters/ nfreq,alpha_s, &
+					beta_s,Tmin,Tmax,gamma_s
+	namelist /model_parameter/ nbin,dx,Cfl
+	namelist /waves_parameters/ Tm,Hs,disp
+	namelist /ice_parameters/ cice,hice,D0,gam,Dmin
+
+	open(30,file='nml/parameter.nml',status='old')
+	read(30,nml=spectrum_parameters)	
+ 	close(30)
+
+	open(30,file='nml/parameter.nml',status='old')
+	read(30,nml=model_parameter)	
+ 	close(30)
+
+	open(30,file='nml/parameter.nml',status='old')
+	read(30,nml=waves_parameters)	
+ 	close(30)
+
+	open(30,file='nml/parameter.nml',status='old')
+	read(30,nml=ice_parameters)	
+ 	close(30)
+	
+
+end subroutine read_namelist
+
+
+
+subroutine array_allocation
+
+nfreq=nfreq+1
+
+
+
+
+allocate(Ei(nfreq))
+allocate(T(nfreq))
+allocate(omega(nfreq))
+allocate(freq(nfreq))
+allocate(sigma_s(nfreq))
+allocate(wl(nfreq))
+allocate(Cp(nfreq))
+allocate(Cg(nfreq))
+allocate(CN(nfreq))
+allocate(alpha(nbin,nfreq))
+allocate(h(nbin))
+allocate(C_ice(nbin))
+allocate(S_ice(nbin,nfreq))
+
+
+	
+	domega=(2*pi/Tmin-2*pi/Tmax)/(nfreq-1)
+	ii=1
+	do i=0,nfreq-1
+	omega(ii)=2*pi/Tmax+i*domega
+	freq(ii)=1/Tmax+i*(1/Tmin-1/Tmax)/(nfreq-1)
+	ii=ii+1
+	end do
+	T=1/freq
+	freq_s=1/Tm
+	wl=g*T**2/(2*pi)
+	Cp=sqrt(g*wl/(2*pi))
+	Cg=Cp/2
+
+	
+	
+	if (disp.eq.0) then
+	CN=Cfl
+	Cg=maxval(Cg)
+	dt=CN(1)*dx/maxval(Cg)	
+	end if
+
+	
+
+	nsteps=ceiling(nbin/minval(CN))
+
+	h(1:nbin)=hice
+	C_ice(1:nbin)=cice
+
+	
+
+	
+allocate(E(nsteps,nbin,nfreq))	
+allocate(Dmax(nbin))
+allocate(Dave(nbin))
+
+	Dmax=D0
+	Dave=D0
+	Dmax(1)=0
+	Dave(1)=0
+
+end subroutine array_allocation
+
+
+	
+
+
+
+
+end module parameters

src/write_output.f90

+subroutine