Skip to content

GitLab

W
WIM2
Commit 37b993a2 authored by Dany Dumont's avatar Dany Dumont
Browse files
Options

Nettoyage et resolution du bogue netcdf

parent 31a3ff4d
Hide whitespace changes
Inline Side-by-side
Showing with 182 additions and 154 deletions
nml/parameter.nml
View file @ 37b993a2
......@@ -20,7 +20,7 @@
Tm =4.5
Hs =1.5
disp =0
disp =1
cont =1
/
......@@ -37,8 +37,8 @@ cont =1
!--------------------------------------------------------------------
&model_parameter
nx =1000
dx =5
nx =100
dx =10
cfl =1.0
name_sim ='test'
root = 'output/'
......@@ -91,7 +91,7 @@ swell_Hs =0.09
X_ice =50
cice =1
ice_thick=1
hice =0
hice =0.8
hmax =0.8
Xh =250
strain_crit=3e-5
......@@ -131,7 +131,7 @@ nr =25
!--------------------------------------------------------------------
&itd_parameters
itd_scheme =1
itd_scheme =0
mu_itd =0.5
mincat_h =0.1
maxcat_h =10
......
src/attenuation.f90
View file @ 37b993a2
!____________________________________________________________________
!
!--------------------------------------------------------------------
!DESCRIPTION: In this routine, the attenuation coefficient is
! calculated according Kohout and Meylan (2008) and the
! spectrum attenuation is 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 (nf))
allocate(p2 (nf))
allocate(p3 (nf))
allocate(att(nf))
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*(p1*h(i)**2 + p2*h(i) + p3)
do ii=1,nf
if ( alpha(i,ii).lt.0 )then
alpha(i,ii)=0d0
end if
end do
att=C_ice(i)*alpha(i,:)/(Dave(i)+3e-14)
!S_ice(n,i,1:nf)=-att*E(n,i,1:nf)
S_ice(i,:,n)=-att*E(i,:,n)
E(i,:,n)=E(i,:,n)*exp(-att*Cg*dt)
!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 (nf))
allocate(p2 (nf))
allocate(p3 (nf))
allocate(att(nf))
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*(p1*h(i)**2 + p2*h(i) + p3)
do ii=1,nf
if ( alpha(i,ii).lt.0 )then
alpha(i,ii)=0d0
end if
end do
att=C_ice(i)*alpha(i,:)/(Dave(i)+3e-14)
!S_ice(n,i,1:nf)=-att*E(n,i,1:nf)
S_ice(i,:,n)=-att*E(i,:,n)
E(i,:,n)=E(i,:,n)*exp(-att*Cg*dt)
end subroutine attenuation
src/ice_fracture.f90
View file @ 37b993a2
!--------------------------------------------------------------------
!DESCRIPTION: This routine does the wave-induced floe breaking by
! modifying the floe size distribution.
!--------------------------------------------------------------------
! There is a major issue with as the integral of the distribution
! is not consistent with the ice area and the redistribution is not
! performed conservatively.
subroutine ice_fracture
......@@ -6,20 +15,18 @@ subroutine ice_fracture
!local parameters
implicit none
integer :: Nband=500
!integer :: Nband=500
integer :: freq1
integer :: freq2
integer :: ninterp=1
double precision :: Y=5.5e9
double precision :: TT=13
double precision :: poisson=0.3
!double precision :: TT=13
double precision :: m_0
double precision :: m_0_stress
double precision :: m_0_amp
double precision :: m_2_amp
double precision :: Tw
double precision :: kw
!double precision :: Tw
!double precision :: kw
double precision :: Lmin
double precision, allocatable :: P(:)
......@@ -70,9 +77,13 @@ subroutine ice_fracture
do iii=1,nh
! Lmin according Mellor 1986
Lmin=pi*0.5*((5e6*middle_h_cat(iii)**3)/(3*10*(1-0.3**2)))**0.25
W=(9.81*kice(iii,:))/omega**2
!Lmin=pi*0.5*((Y*middle_h_cat(iii)**3)/ &
! (3*rhoi*g*(1-nu**2)))**0.25
! Lmin corrected from Mellor (1986)
Lmin=pi*0.25*((Y*middle_h_cat(iii)**3)/ &
(3*rhoi*g*(1-nu**2)))**0.25
W=(g*kice(iii,:))/omega**2
wavelength=2*pi/kice(iii,:)
strain=middle_h_cat(iii)*0.5*kice(iii,:)**2*W
......@@ -82,8 +93,8 @@ subroutine ice_fracture
do ii=1,nr
if ( 2d0*floe_cat(ii+1).lt.wavelength(nf) .or. &
2d0*floe_cat(ii+1).gt.wavelength(1)) then
if ( 2d0*floe_cat(ii+1).lt.wavelength(nf) .or. &
2d0*floe_cat(ii+1).gt.wavelength( 1) ) then
P_cat(ii)=0d0
else
do while ( wavelength(freq1).lt.2d0*floe_cat(ii+1) .or. &
......
src/initialization.f90
View file @ 37b993a2
......@@ -42,6 +42,7 @@ subroutine initialization
end do
write(*,*) ' * dt = ',dt,'s'
write(*,*) ' * dx = ',dx,'m'
if ( disp.eq.1 ) then
write(*,*) ' * dispersion is on'
else
......@@ -49,7 +50,7 @@ subroutine initialization
end if
if ( cont.eq.1 ) then
write(*,*) ' * wave energy is continuously injected'
write(*,*) ' * wave energy continuously injected'
else
write(*,*) ' * only a single wave energy pulse'
end if
......@@ -62,6 +63,8 @@ subroutine initialization
if ( init_spec.eq.1 ) then
write(*,*) ' * Jonswap spectrum'
! JONSWAP spectrum
do i=1,nf
if ( freq(i).le.freq_s ) then
......@@ -72,18 +75,22 @@ subroutine initialization
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)
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)* &
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=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
......@@ -94,12 +101,12 @@ subroutine initialization
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( 1:x1)=0 !ice concentration is 0 before ice edge
C_ice(x1:nx)=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:nx)=maxfloe !initalize mean floe size in ice 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
......@@ -118,10 +125,10 @@ subroutine initialization
!--------------------- FLOE SIZE DISTRIBUTION -----------------------
res=abs(minfloe-maxfloe)/nedge
dr=abs(minfloe - maxfloe)/nedge
floe_cat(1)=minfloe
do i=2,nedge
floe_cat(i)=floe_cat(1)+i*res
floe_cat(i)=floe_cat(1)+i*dr
end do
do i=1,nr
......@@ -134,18 +141,17 @@ subroutine initialization
end do
Dave(x1:nx)=middle_floe_cat(nr)
!initialize max floe size in ice transect
!--------------------- ICE THICKNESS DISTRIBUTION -------------------
resh=abs(mincat_h-maxcat_h)/(nedge_h)
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*resh
h_cat(i)=h_cat(1) + i*dh
end do
do i=1,nh
middle_h_cat(i)=(h_cat(i)+h_cat(i+1))*0.5
middle_h_cat(i)=0.5*(h_cat(i) + h_cat(i+1))
end do
itd(:,:)=0d0
......@@ -170,7 +176,6 @@ subroutine initialization
end if
!--------------------- STATISTICS -----------------------------------
h_sign(1,1)=Hs
......
src/main.f90
View file @ 37b993a2
!
!____________________________________________________________________
!
!--------------------------------------------------------------------
!DESCRIPTION: This is the main program of WIM. This routine merely
! calls other subroutines and do the main time loop.
! It also contains the subroutine progress which display
! a progress bar in the terminal while the model is
! running.
!____________________________________________________________________
!--------------------------------------------------------------------
!INTERFACE:
program wim2
......@@ -113,7 +111,7 @@ subroutine message(info)
endif
if ( info.eq.3 ) then
write(*,*)' Allocate memory... '
write(*,*)' Allocate memory ... '
endif
if ( info.eq.4 ) then
......@@ -129,12 +127,9 @@ subroutine message(info)
endif
if ( info.eq.7 ) then
write(*,*)' done!'
write(*,*)' done!'
endif
end subroutine message
!____________________________________________________________________
end program wim2
src/parameters.f90
View file @ 37b993a2
!____________________________________________________________________
!
!--------------------------------------------------------------------
!DESCRIPTION: In this routine, the shared parameters for other
! subroutines are declared. It also contains the
! subroutine read_namelist which read parameters value
! from the namelist and the subroutine array_allocation
! for memory allocation.
!____________________________________________________________________
!--------------------------------------------------------------------
!INTERFACE:
module parameters
implicit none
!NAME OF OUTPUTS FILES
character(len=100) :: root='output/' ! Output directory
character(len=100) :: name_sim='test'! Name of the simulation
character(len=100) :: namefile
character(len=100) :: dirout
character(len=100) :: root='output/' ! Output directory
character(len=100) :: name_sim='test'! Name of the simulation
character(len=100) :: namefile
character(len=100) :: dirout
!DUMMIES
integer :: i,ii,iii
integer :: j,jj,jjj
integer :: n
integer :: i,ii,iii
integer :: j,jj,jjj
integer :: n
!GLOBAL PARAMETERS
double precision :: g=9.81 ! Gravitationnal acceleration
double precision :: pi=3.1416 ! Pi
! physical and useful quantities
real :: g=9.81 ! Gravitationnal acc. [m s-2]
real :: Y=5.5e9 ! Young's modulus [Pa]
real :: nu=0.3 ! Poisson's ratio [-]
real :: rhoi=925.0 ! Sea ice density [kg m-3]
double precision :: pi=3.141592653 ! Pi
!SPECTRUM
integer :: nf=60 ! Nb of frequency bins
integer :: init_spec=1
double precision :: alpha_s=8.1e-3
double precision :: beta_s=1.25
double precision :: Tmin=2.5 ! Minimum period
double precision :: Tmax=20 ! Maximum period
double precision :: Tm=6 ! Peak period
double precision :: Hs=1 ! Significant Height
double precision :: gamma_s=3.3
double precision :: freq_s ! Peak frequency
double precision :: domega ! Frequency increment
double precision :: swell_T=8
double precision :: swell_Hs=1
integer :: nf=60 ! Nb of frequency bins
integer :: init_spec=1
double precision :: alpha_s=8.1e-3 ! Jonswap parameter
double precision :: beta_s=1.25 ! Jonswap parameter
double precision :: gamma_s=3.3 ! Jonswap parameter
double precision :: Tmin=2.5 ! Minimum period
double precision :: Tmax=20 ! Maximum period
double precision :: Tm=6 ! Peak period [s]
double precision :: Hs=1 ! Significant wave height [m]
double precision :: freq_s ! Peak frequency
double precision :: domega ! Frequency increment
double precision :: swell_T=8
double precision :: swell_Hs=1
double precision, allocatable :: E(:,:,:) ! Energy spectrum
double precision, allocatable :: Ei(:) ! Initial spectrum
......@@ -67,7 +67,7 @@ integer :: nx=10 ! Nb of spatial bins
integer :: nt ! Duration of the simulation
real :: dt ! Time step
double precision :: dx=500 ! Grid resolution
real :: dx=500 ! Grid resolution
double precision :: cfl=1 ! Courant
double precision, allocatable :: x_axis(:) ! Vector of grid width
......@@ -106,7 +106,7 @@ integer :: nedge
double precision :: minfloe=5
double precision :: maxfloe=500
double precision :: res
double precision :: dr
double precision, allocatable :: fsd(:,:,:,:)
double precision, allocatable :: floe_cat(:)
......@@ -119,8 +119,8 @@ integer :: nedge_h
double precision :: mu_itd=1.2
double precision :: mincat_h=0.1
double precision :: maxcat_h=10
double precision :: resh
double precision :: maxcat_h=5.0
double precision :: dh
double precision, allocatable :: itd(:,:)
double precision, allocatable :: middle_h_cat(:)
......@@ -135,7 +135,7 @@ contains
subroutine read_namelist
namelist /spectrum_parameters/ nf,alpha_s,beta_s, &
namelist /spectrum_parameters/ nf,alpha_s,beta_s, &
Tmin,Tmax,gamma_s,init_spec, &
swell_T,swell_Hs
......@@ -195,7 +195,7 @@ subroutine array_allocation
allocate(C_ice (nx))
allocate(x_axis (nx))
domega=(2*pi/Tmin-2*pi/Tmax)/(nf)
domega=(2*pi/Tmin - 2*pi/Tmax)/nf
dwl=(g*Tmax**2/(2*pi)-g*Tmin**2/(2*pi))/(nf)
ii=1
do i=0,nf
......@@ -222,11 +222,11 @@ subroutine array_allocation
nt=ceiling(nx/minval(CN))
allocate(E (nx,nf,nt))
allocate(Dmax (nx))
allocate(Dave (nx))
allocate(time (nt))
allocate(Dmax (nx ))
allocate(Dave (nx ))
allocate(time (nt ))
allocate(S_ice (nx,nf,nt))
allocate(h_sign(nx,nt))
allocate(h_sign(nx,nt ))
nedge=nr+1
......
src/write_output.f90
View file @ 37b993a2
......@@ -22,11 +22,13 @@ subroutine write_output
include 'netcdf.inc'
integer :: ncid
integer :: omID
integer :: xID
integer :: fID
integer :: rID
integer :: hID
integer :: tID
integer :: EID
integer :: om_varID
integer :: f_varID
integer :: x_varID
integer :: t_varID
integer :: stat
......@@ -34,12 +36,10 @@ subroutine write_output
integer :: Dave_varID
integer :: c_varID
integer :: h_varID
integer :: floeID
integer :: fsdID
integer :: floe_varID
integer :: r_varID
integer :: SrcID
integer :: HsID
integer :: hID
integer :: hcat_varID
integer :: itdID
......@@ -47,7 +47,7 @@ subroutine write_output
call handle_err(stat)
! dimensions
stat=nf90_def_dim(ncid,"omega",nf,omID)
stat=nf90_def_dim(ncid,"omega",nf,fID)
call handle_err(stat)
stat=nf90_def_dim(ncid,"x" ,nx,xID)
......@@ -56,14 +56,14 @@ subroutine write_output
stat=nf90_def_dim(ncid,"time",nt,tID)
call handle_err(stat)
stat=nf90_def_dim(ncid,"dcat",nr,floeID)
stat=nf90_def_dim(ncid,"dcat",nr,rID)
call handle_err(stat)
stat=nf90_def_dim(ncid,"hcat",nh,hID)
call handle_err(stat)
! coordinates
stat=nf90_def_var(ncid,"omega",nf90_double,omID,om_varID)
stat=nf90_def_var(ncid,"omega",nf90_double,fID,f_varID)
call handle_err(stat)
stat=nf90_def_var(ncid,"x",nf90_double,xID,x_varID)
......@@ -75,15 +75,15 @@ subroutine write_output
stat=nf90_def_var(ncid,"hcat",nf90_double,hID,hcat_varID)
call handle_err(stat)
stat=nf90_def_var(ncid,"dcat",nf90_double,floeID,floe_varID)
stat=nf90_def_var(ncid,"dcat",nf90_double,rID,r_varID)
call handle_err(stat)
! variables
! (x, f, t)
stat=nf90_def_var(ncid,"E",nf90_double,(/xID,omID,tID/),EID)
stat=nf90_def_var(ncid,"E",nf90_double,(/xID,fID,tID/),EID)
call handle_err(stat)
stat=nf90_def_var(ncid,"Sice",nf90_double,(/xID,omID,tID/),SrcID)
stat=nf90_def_var(ncid,"Sice",nf90_double,(/xID,fID,tID/),SrcID)
call handle_err(stat)
! (x)
......@@ -100,7 +100,7 @@ subroutine write_output
call handle_err(stat)
! (x, r, h, t)
stat=nf90_def_var(ncid,"fsd",nf90_double,(/xID,floeID,hID,tID/),fsdID)
stat=nf90_def_var(ncid,"fsd",nf90_double,(/xID,rID,hID,tID/),fsdID)
call handle_err(stat)
! (x, t)
......@@ -113,37 +113,58 @@ subroutine write_output
! set attributes
stat=nf_put_att_text(ncid,xID,'units',10,'km')
call handle_err(stat)
stat=nf_put_att_text(ncid,xID,'axis',1,'X')
stat=nf_put_att_text(ncid,xID,'long_name',7,'Distance')
call handle_err(stat)
stat=nf_put_att_text(ncid,xID,'long_name',10,'Distance')
call handle_err(stat)
stat=nf_put_att_text(ncid,omID,'units',5,'Hz')
stat=nf_put_att_text(ncid,omID,'long_name',7,'Frequency')
stat=nf_put_att_text(ncid,rID,'units',5,'m')
call handle_err(stat)
stat=nf_put_att_text(ncid,rID,'long_name',10,'Floe size')
call handle_err(stat)
stat=nf_put_att_text(ncid,fID,'units',5,'Hz')
call handle_err(stat)
stat=nf_put_att_text(ncid,fID,'long_name',10,'Frequency')
call handle_err(stat)
stat=nf_put_att_text(ncid,hID,'units',5,'m')
call handle_err(stat)
stat=nf_put_att_text(ncid,hID,'axis',1,'Z')