Commit 31a3ff4d authored by Dany Dumont's avatar Dany Dumont

Changement de la structure de la variable fsd (t,x,r,h) > (x,r,h,t)

parent c2c1e911
......@@ -20,7 +20,7 @@
Tm =4.5
Hs =1.5
disp =1
disp =0
cont =1
/
......@@ -29,7 +29,7 @@ cont =1
!
! nx -> Number of grid bin
! dx -> Spatial resolution [m]
! Cfl -> Courant-Friedrich-Lewy condition (0 < Cfl < 1)
! cfl -> Courant-Friedrich-Lewy condition (0 < cfl < 1)
! Only in the case where disp=0. The CFL condition
! is needed to calculate the time step.
! name_sim -> name of the output file
......@@ -37,9 +37,9 @@ cont =1
!--------------------------------------------------------------------
&model_parameter
nx =100
nx =1000
dx =5
Cfl =1.0
cfl =1.0
name_sim ='test'
root = 'output/'
/
......@@ -62,7 +62,7 @@ root = 'output/'
!--------------------------------------------------------------------
&spectrum_parameters
init_spec =0
init_spec =1
nf =36
Tmin =2.5
Tmax =20
......@@ -88,9 +88,9 @@ swell_Hs =0.09
!--------------------------------------------------------------------
&ice_parameters
X_ice =10
X_ice =50
cice =1
ice_thick=0
ice_thick=1
hice =0
hmax =0.8
Xh =250
......@@ -101,7 +101,7 @@ P_c =0.55
!____________________________________________________________________
! FSD PARAMETERS
! FSD_sheme -> method for compute <D>
! fsd_sheme -> method for compute <D>
! 0: dumont et al (2011)
! 1: power law
!
......@@ -111,10 +111,10 @@ P_c =0.55
!--------------------------------------------------------------------
&fsd_parameters
FSD_scheme =1
fsd_scheme =1
minfloe =5
maxfloe =400
nr =60
nr =25
/
!____________________________________________________________________
......
......@@ -129,11 +129,14 @@ subroutine ice_fracture
!compute the new floe size distribution!
do jj=1,nr
Q2(jj)=sum(Q(jj,:)*FSD(n-1,i,:,iii)*F)
!Q2(jj)=sum(Q(jj,:)*fsd(n-1,i,:,iii)*F)
Q2(jj)=sum(Q(jj,:)*fsd(i,:,iii,n-1)*F)
end do
FSD(n,i,1:nr,iii)=FSD(n-1,i,1:nr,iii) - &
F*FSD(n-1,i,1:nr,iii) + Q2
!fsd(n,i,1:nr,iii)=fsd(n-1,i,1:nr,iii) - &
! F*fsd(n-1,i,1:nr,iii) + Q2
fsd(i,1:nr,iii,n)=fsd(i,1:nr,iii,n-1) - &
F*fsd(i,1:nr,iii,n-1) + Q2
!if (i.eq.25 .and.n.eq.25.and.iii.eq.10)then
!open(31,file='output/strain3.dat')
......@@ -150,8 +153,8 @@ subroutine ice_fracture
!compute the new average floe diameter <D>
do ii=1,nh
nfloe(:,ii)=FSD(n,i,:,ii)*itd(i,ii)*dx**2/middle_floe_cat**2
Dave1(ii)=sum(FSD(n,i,:,ii)*middle_floe_cat)
nfloe(:,ii)=fsd(i,:,ii,n)*itd(i,ii)*dx**2/middle_floe_cat**2
Dave1(ii)=sum(fsd(i,:,ii,n)*middle_floe_cat)
enddo
nfloe=nfloe/sum(nfloe+3e-14)
......
......@@ -28,20 +28,33 @@ subroutine initialization
allocate(Gf(nf))
allocate(PM(nf))
! construct time array in hours
! construct time array [in hours]
time(1)=0
do ii=2,nt
time(ii)=time(ii-1)+dt/3600
end do
write(*,*) ' * Time step is ',dt,' s'
! construct space array in km
! construct space array [in km]
x_axis(1)=0
do ii=2,nx
x_axis(ii)=x_axis(ii-1)+dx/1000
x_axis(ii)=x_axis(ii-1) + dx/1000
end do
write(*,*) ' * dt = ',dt,'s'
if ( disp.eq.1 ) then
write(*,*) ' * dispersion is on'
else
write(*,*) ' * dispersion is off'
end if
if ( cont.eq.1 ) then
write(*,*) ' * wave energy is continuously injected'
else
write(*,*) ' * only a single wave energy pulse'
end if
!--------------------- INITIAL SPECTRUM -----------------------------
......@@ -75,11 +88,11 @@ subroutine initialization
end if
E(1,1:nf,1)=Ei
E(1,:,1)=Ei
!--------------------- ICE TRANSECT ---------------------------------
S_ice(1,1:nf,1)=0
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
......@@ -91,7 +104,7 @@ subroutine initialization
if ( ice_thick.eq.0 ) then
! constant ice thickness in the transect
h(x1:nx)=hice
h(x1:nx )=hice
h(1 :x1-1)=0d0
else
......@@ -115,9 +128,9 @@ subroutine initialization
middle_floe_cat(i)=(floe_cat(i)+floe_cat(i+1))*0.5
end do
FSD(:,:,:,:)=0d0
fsd(:,:,:,:)=0d0
do i=x1,nx
FSD(1,i,nr,:)=C_ice(i)
fsd(i,nr,:,1)=C_ice(i)
end do
Dave(x1:nx)=middle_floe_cat(nr)
......
......@@ -43,7 +43,7 @@ call message(5)
do n=2,nt
! display progress bar
call progress(n,nt)
!call progress(n,nt)
! do advection for each frequency band
do ii=1,nf
......@@ -79,9 +79,9 @@ subroutine progress(j,jmax)
implicit none
integer :: j,k,jmax
character(len=39) :: bar=' processing ???%'
character(len=20) :: bar=' processing ???%'
write(unit=bar(36:38),fmt="(i3)")ceiling((real(j)/real(jmax))*100)
write(unit=bar(17:19),fmt="(i3)")ceiling((real(j)/real(jmax))*100)
write(unit=6,fmt="(a1,a22)",advance="no") char(13), bar
if (real(j)/real(jmax).eq.1) then
......
......@@ -68,7 +68,7 @@ integer :: nt ! Duration of the simulation
real :: dt ! Time step
double precision :: dx=500 ! Grid resolution
double precision :: Cfl=1 ! Courant
double precision :: cfl=1 ! Courant
double precision, allocatable :: x_axis(:) ! Vector of grid width
double precision, allocatable :: time(:) ! Vector of time
......@@ -100,7 +100,7 @@ double precision, allocatable :: Dmax(:) ! Maximum floe size
double precision, allocatable :: kice(:,:)
!FSD
integer :: FSD_scheme=1
integer :: fsd_scheme=1
integer :: nr=40
integer :: nedge
......@@ -108,7 +108,7 @@ double precision :: minfloe=5
double precision :: maxfloe=500
double precision :: res
double precision, allocatable :: FSD(:,:,:,:)
double precision, allocatable :: fsd(:,:,:,:)
double precision, allocatable :: floe_cat(:)
double precision, allocatable :: middle_floe_cat(:)
......@@ -139,14 +139,14 @@ subroutine read_namelist
Tmin,Tmax,gamma_s,init_spec, &
swell_T,swell_Hs
namelist /model_parameter/ nx,dx,Cfl,name_sim,root
namelist /model_parameter/ nx,dx,cfl,name_sim,root
namelist /waves_parameters/ Tm,Hs,disp,cont
namelist /ice_parameters/ cice,hice,ice_thick,hmax, &
X_ice,Xh,strain_crit,P_c
namelist /fsd_parameters/ FSD_scheme,minfloe,maxfloe, &
namelist /fsd_parameters/ fsd_scheme,minfloe,maxfloe, &
nr
namelist /itd_parameters/ itd_scheme,mu_itd,mincat_h, &
......@@ -211,7 +211,7 @@ subroutine array_allocation
Cg=0.5*Cp
if ( disp.eq.0 ) then
CN=Cfl
CN=cfl
Cg=maxval(Cg)
dt=CN(1)*dx/maxval(Cg)
else
......@@ -226,12 +226,12 @@ subroutine array_allocation
allocate(Dave (nx))
allocate(time (nt))
allocate(S_ice (nx,nf,nt))
allocate(h_sign(nt,nx))
allocate(h_sign(nx,nt))
nedge=nr+1
allocate(floe_cat (nedge))
allocate(FSD (nt,nx,nr,nh))
allocate(fsd (nx,nr,nh,nt))
allocate(middle_floe_cat(nr))
nedge_h=nh+1
......
......@@ -7,6 +7,6 @@ subroutine statistics
double precision :: m_0
m_0=sum(E(i,:,n)*domega)
h_sign(n,i)=4*sqrt(m_0)
h_sign(i,n)=4*sqrt(m_0)
end subroutine
......@@ -41,7 +41,7 @@ subroutine write_output
integer :: HsID
integer :: hID
integer :: hcat_varID
integer :: idtID
integer :: itdID
stat=nf90_create(namefile,cmode=or(nf90_clobber,nf90_64bit_offset),ncid=ncid)
call handle_err(stat)
......@@ -99,16 +99,16 @@ subroutine write_output
stat=nf90_def_var(ncid,"ithk",nf90_double,xID,h_varID)
call handle_err(stat)
! (t, x, d, h)
stat=nf90_def_var(ncid,"fsd",nf90_double,(/tID,xID,floeID,hID/),fsdID)
! (x, r, h, t)
stat=nf90_def_var(ncid,"fsd",nf90_double,(/xID,floeID,hID,tID/),fsdID)
call handle_err(stat)
! (t, x)
stat=nf90_def_var(ncid,"hs",nf90_double,(/tID,xID/),HsID)
! (x, t)
stat=nf90_def_var(ncid,"hs",nf90_double,(/xID,tID/),HsID)
call handle_err(stat)
! (x, h)
stat=nf90_def_var(ncid,"itd",nf90_double,(/xID,hID/),idtID)
stat=nf90_def_var(ncid,"itd",nf90_double,(/xID,hID/),itdID)
call handle_err(stat)
! set attributes
......@@ -135,6 +135,9 @@ subroutine write_output
stat=nf_put_att_text(ncid,EID,'long_name',50, &
'Sea surface wave variance spectral density')
stat=nf_put_att_text(ncid,fsdID,'units',10,'m m-1')
stat=nf_put_att_text(ncid,fsdID,'long_name',50, &
'Floe size and thickness distribution')
! end file definition
stat=nf90_enddef(ncid)
call handle_err(stat)
......@@ -167,7 +170,7 @@ subroutine write_output
stat=nf90_put_var(ncid,h_varID,h)
call handle_err(stat)
stat=nf90_put_var(ncid,fsdID,FSD)
stat=nf90_put_var(ncid,fsdID,fsd)
call handle_err(stat)
stat=nf90_put_var(ncid,floe_varID,middle_floe_cat)
......@@ -179,7 +182,7 @@ subroutine write_output
stat=nf90_put_var(ncid,HsID,h_sign)
call handle_err(stat)
stat=nf90_put_var(ncid,idtID,itd)
stat=nf90_put_var(ncid,itdID,itd)
call handle_err(stat)
stat=nf90_close(ncid)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment