Commit 4e5b4414 authored by Jérémy Baudry's avatar Jérémy Baudry
Browse files

ajout commentaires dans le code

parent a081b91c
...@@ -17,7 +17,7 @@ ...@@ -17,7 +17,7 @@
&waves_parameters &waves_parameters
Tm =6 Tm =6
Hs =1.5 Hs =1
disp =0 disp =0
/ /
...@@ -31,14 +31,14 @@ disp =0 ...@@ -31,14 +31,14 @@ disp =0
! is needed to calculate the time step. ! is needed to calculate the time step.
! name_sim -> name of the output file ! name_sim -> name of the output file
! root -> destination folder for the output file ! root -> destination folder for the output file
! FSD_sheme -> method for ice floe breaking parametrization ! FSD_sheme -> method for compute <D>
! 0: Williams et al. 2013 ! 0: dumont et al (2011)
! 1: Horvat 2015 ! 1: power law
!------------------------------------------------------------------------------ !------------------------------------------------------------------------------
&model_parameter &model_parameter
nbin =100 nbin =50
dx =1000 dx =5000
Cfl =1 Cfl =1
name_sim ='simulation1' name_sim ='simulation1'
root = 'output/' root = 'output/'
...@@ -86,9 +86,9 @@ gamma_s =3.3 ...@@ -86,9 +86,9 @@ gamma_s =3.3
! Dmin -> Minimum floe size (if FSD_sheme=1) [m] ! Dmin -> Minimum floe size (if FSD_sheme=1) [m]
!------------------------------------------------------------------------------ !------------------------------------------------------------------------------
&ice_parameters &ice_parameters
X_ice =5000 X_ice =50000
cice =0.8 cice =0.8
ice_thick =0 ice_thick =1
hice =2 hice =2
hmax =4 hmax =4
Xh =100000 Xh =100000
......
subroutine advection
use parameters !_______________________________________________________________________________
implicit none !DESCRIPTION: Here, the wave spectrum is advected through the
!domain. The advection equation is solved using a Lax-wendroff
!discretization sheme with a superbee flux limiter.
double precision, allocatable :: diffF(:),diffl(:),diffr(:),phi(:),theta(:) !_______________________________________________________________________________
double precision, allocatable :: F(:),diff1(:)
!INTERFACE:
subroutine advection
!MODULE USES:
use parameters
!LOCAL PARAMETERS:
implicit none
double precision, allocatable :: diffF(:)
double precision, allocatable :: diffl(:)
double precision, allocatable :: diffr(:)
double precision, allocatable :: phi(:)
double precision, allocatable :: theta(:)
double precision, allocatable :: F(:)
double precision, allocatable :: diff1(:)
allocate(diffl(nbin)) allocate(diffl(nbin))
allocate(diffr(nbin)) allocate(diffr(nbin))
allocate(diffF(nbin)) allocate(diffF(nbin))
allocate(theta(nbin)) allocate(theta(nbin))
allocate(phi(nbin)) allocate(phi(nbin))
allocate(F(nbin)) allocate(F(nbin))
allocate(diff1(nbin-1)) allocate(diff1(nbin-1))
!________________________________________________________________________________
do i=2,nbin do i=2,nbin
diffl(i)=E(n-1,i,ii)-E(n-1,i-1,ii) diffl(i)=E(n-1,i,ii)-E(n-1,i-1,ii)
end do end do
do i=1,nbin-1 do i=1,nbin-1
diffr(i)=E(n-1,i+1,ii)-E(n-1,i,ii) diffr(i)=E(n-1,i+1,ii)-E(n-1,i,ii)
end do end do
diffr(nbin)=-E(n-1,nbin,ii)
diffr(nbin)=-E(n-1,nbin,ii) !set dE=0 at boundaries
diffl(1)=E(n-1,1,ii) diffl(1)=E(n-1,1,ii)
do i=1,nbin do i=1,nbin !Superbee flux limiter
theta(i)=diffl(i)/(diffr(i)+3e-14) theta(i)=diffl(i)/(diffr(i)+3e-14)
phi(i)=max(0d0,min(theta(i)*2,1d0),min(theta(i),2d0)) phi(i)=max(0d0,min(theta(i)*2,1d0),min(theta(i),2d0))
end do end do
!Lax-Wendroff sheme:
F=E(n-1,1:nbin,ii)+0.5*(1-CN(ii))*diffr*phi F=E(n-1,1:nbin,ii)+0.5*(1-CN(ii))*diffr*phi
do i=2,nbin do i=2,nbin
diffF(i)=F(i)-F(i-1) diffF(i)=F(i)-F(i-1)
...@@ -42,7 +63,7 @@ implicit none ...@@ -42,7 +63,7 @@ implicit none
E(n,1:nbin,ii)=E(n-1,1:nbin,ii)-CN(ii)*diffF E(n,1:nbin,ii)=E(n-1,1:nbin,ii)-CN(ii)*diffF
!________________________________________________________________________________
end subroutine advection end subroutine advection
......
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)) !DESCRIPTION: In this routine, the attenuation
allocate(p2(nfreq)) !coefficient is calculated according Kohout and
allocate(p3(nfreq)) !Meylan (2008) and the spectrum attenuation is
allocate(att(nfreq)) !computed.
!____________________________________________________________________________
!INTERFACE:
subroutine attenuation
!USES:
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(:)
double precision, allocatable :: att(:)
allocate(p1(nfreq))
allocate(p2(nfreq))
allocate(p3(nfreq))
allocate(att(nfreq))
!____________________________________________________________________________
q11 = -0.00000777 q11 = -0.00000777
q12 = 0.00032080 q12 = 0.00032080
...@@ -42,13 +60,14 @@ alpha(i,j)=0d0 ...@@ -42,13 +60,14 @@ alpha(i,j)=0d0
end if end if
end do end do
att=C_ice(i)*alpha(i,1:nfreq)/(Dave(i)+3e-14) att=C_ice(i)*alpha(i,1:nfreq)/(Dave(i)+3e-14)
S_ice(i,1:nfreq)=-att*E(n,i,1:nfreq) S_ice(n,i,1:nfreq)=-att*E(n,i,1:nfreq)
E(n,i,1:nfreq)=E(n,i,1:nfreq)*exp(-att*Cg*dt) E(n,i,1:nfreq)=E(n,i,1:nfreq)*exp(-att*Cg*dt)
!_____________________________________________________________________________
end subroutine attenuation end subroutine attenuation
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=10000
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/12d0*(1-v**2)
d=0.9*h(i)
jj=1
do ii=0,acc
k(jj)=real(ii)*10d0/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=(0.5*h(i))*WN_ice**2*W !DESCRIPTION: Here is the parametrization proposed by
m0_s=sum(E(n,i,1:nfreq)*S**2)*domega !Williams et al. (2013) for wave induced
SS=2*sqrt(m0_s) !floe breaking. In this routine, the maximum floe size
if (SS.gt.RS) then !is computed using the significant strain imposed by
!passing waves in ice.
!_______________________________________________________________________________
!INTERFACE:
subroutine floe_breaking
!USES:
use parameters
!LOCAL PARAMETERS:
implicit none
m0_a=sum(E(n,i,1:nfreq)*W**2)*domega double precision, allocatable :: WN_ice(:)
m2_a=sum(omega**2*E(n,i,1:nfreq)*W**2)*domega double precision, allocatable :: fnt(:)
Tw=2*pi*sqrt(m0_a/m2_a) double precision, allocatable :: k(:)
om_d=2*pi/Tw double precision :: rho=1025
double precision :: d
double precision :: F
double precision :: Y=5.5
double precision :: v=0.3
integer :: acc=10000
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/12d0*(1-v**2)
d=0.9*h(i)
jj=1 !create a wavelength vector (this is used to
!find the root of the polynom fnt
do ii=0,acc
k(jj)=real(ii)*10d0/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))!calculate wavelength in ice
end do
W=g*WN_ice/omega**2
S=(0.5*h(i))*WN_ice**2*W
m0_s=sum(E(n,i,1:nfreq)*S**2)*domega !0th moment of the strain spectrum
SS=2*sqrt(m0_s) !significant strain
if (SS.gt.RS) then
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 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
Dmax(i)=max(Dmin,min(0.5*wl_w,Dmax(i))) 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(0.5*wl_w,Dmax(i)))
end if
end if
!_________________________________________________________________________________
end subroutine floe_breaking end subroutine floe_breaking
subroutine fsd_build
use parameters
implicit none !________________________________________________________________________________
double precision :: coeff
double precision, allocatable :: ND(:),NN(:) !DESCRIPTION: In this routine the average floe size
integer :: M,mm !Dave is computed using the value of Dmax calculated in
!the subroutine 'floe_breaking'.
!________________________________________________________________________________
!INTERFACE:
subroutine fsd_build
!MODULE USES:
use parameters
!LOCAL PARAMETERS:
implicit none
double precision :: coeff
double precision, allocatable :: ND(:)
double precision, allocatable :: NN(:)
integer :: M
integer :: mm
!_________________________________________________________________________________
if(C_ice(i).eq.0)then if(C_ice(i).eq.0)then
Dave(i)=0 Dave(i)=0
...@@ -13,11 +33,12 @@ use parameters ...@@ -13,11 +33,12 @@ use parameters
elseif (Dmax(i).eq.Dmin) then elseif (Dmax(i).eq.Dmin) then
Dave(i)=Dmin Dave(i)=Dmin
elseif (FSD_scheme.eq.1) then elseif (FSD_scheme.eq.1) then !use a power law to compute <D>
coeff=1/((1/(1-gam))*(Dmax(i)**(1-gam)-Dmin**(1-gam))) 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)) Dave(i)=coeff*(1/(2-gam))*(Dmax(i)**(2-gam)-Dmin**(2-gam))
else
else !use the method from Dumont et al.(2011) to compute <D>
M=floor(log(Dmax(i)/Dmin)/log(psi)) M=floor(log(Dmax(i)/Dmin)/log(psi))
allocate(ND(M+1)) allocate(ND(M+1))
...@@ -40,6 +61,6 @@ use parameters ...@@ -40,6 +61,6 @@ use parameters
!_________________________________________________________________________________
end subroutine fsd_build end subroutine fsd_build
...@@ -23,6 +23,22 @@ ...@@ -23,6 +23,22 @@
!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_____________________________ !_________________________INITIAL SPECTRUM_____________________________
E(1:nsteps,1:nbin,1:nfreq)=0d0 !INITIALIZE SPECTRUM ARRAY E(1:nsteps,1:nbin,1:nfreq)=0d0 !INITIALIZE SPECTRUM ARRAY
...@@ -54,7 +70,7 @@ ...@@ -54,7 +70,7 @@
!_______________________________________________________________________ !_______________________________________________________________________
!_______________________ICE_TRANSECT____________________________________ !_______________________ICE_TRANSECT____________________________________
S_ice(1,1,1:nfreq)=0
X1=floor(X_ice/dx) !find in which grid bin is the ice edge 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(1:X1)=0 !ice concentration is 0 before ice edge!
C_ice(X1:nbin)=cice !ice concentration in the transect C_ice(X1:nbin)=cice !ice concentration in the transect
...@@ -81,6 +97,17 @@ h(1:X1)=0d0 ...@@ -81,6 +97,17 @@ h(1:X1)=0d0
end if end if
!_____________________FSD_____________________________________________ !_____________________FSD_____________________________________________
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
FSD(1,1:nbin,1)=1d0 FSD(1,1:nbin,1)=1d0
......
...@@ -19,18 +19,6 @@ ...@@ -19,18 +19,6 @@
call initialization ! initialize the model call initialization ! initialize the model
!________________________________________________________________________________
spectrum=trim(root)//'Energy_spectrum.dat'
floe_size=trim(root)//'floe_size.dat'
namefile=trim(root)//trim(name_sim)//'.nc'
open(10,file=spectrum)
open(11,file=floe_size)
!________________________________________________________________________________ !________________________________________________________________________________
! DO THE TIME LOOP ! DO THE TIME LOOP
...@@ -46,18 +34,17 @@ ...@@ -46,18 +34,17 @@
do i=1,nbin !spatial calculations do i=1,nbin !spatial calculations
call attenuation ! compute spectrum attenuation call attenuation ! compute spectrum attenuation
call break_horvat ! compute floe breaking ! compute floe breaking
!call floe_breaking call floe_breaking
!call fsd_build call fsd_build
end do end do
end do end do
!______________________OUTPUTS_________________________________________________ !______________________OUTPUTS_________________________________________________
namefile=trim(root)//trim(name_sim)//'.nc'
call write_output ! Write outputs in NETCDF call write_output ! Write outputs in NETCDF
!call graph_dislin
close(10)
close(11)
contains contains
......
...@@ -5,7 +5,7 @@ OPTION= -O3 ...@@ -5,7 +5,7 @@ OPTION= -O3
NETCDFinc=-I/usr/include NETCDFinc=-I/usr/include
NETCDFLIB=-L/usr/lib64 -lnetcdff NETCDFLIB=-L/usr/lib64 -lnetcdff
NETCDFMOD=-I/usr/lib64/gfortran/modules NETCDFMOD=-I/usr/lib64/gfortran/modules
DISLIN=-I/usr/local/dislin/gf -ldislin
OBJ= *.o OBJ= *.o
SRC= *.f90 SRC= *.f90
...@@ -17,7 +17,7 @@ EXEC= WIM2 ...@@ -17,7 +17,7 @@ EXEC= WIM2