Commit b85b16f5 authored by Jérémy Baudry's avatar Jérémy Baudry

clean for winter school

parent d6594dfd
...@@ -11,8 +11,20 @@ ...@@ -11,8 +11,20 @@
! q0 : internal nutrient subsistence quota (for bio_model=2) [mol N/mol C] ! q0 : internal nutrient subsistence quota (for bio_model=2) [mol N/mol C]
! mortality : mortality rate [s-1] ! mortality : mortality rate [s-1]
! Kn : Half saturation constant [mol N/m3] ! Kn : Half saturation constant [mol N/m3]
! Vmax : maximum nutrient uptake rate [molP/molC/s-1] ! Vmax : maximum nutrient uptake rate [molP/molC/s]
! Nutrient_0 : Initial Nutrient concentration [mol N/m3] !
! nut_profile : initial nutrient profile
! 1: uniform in the watercolumn (mixed conditions)
! 2: two layers (stratified conditions)
!
! nut_layer : Depth of the bottom layer (for nut_profile=2) [m]
! Nutrient_0 : Initial Nutrient concentration (for nut_profile=1) [mol N/m3]
! Nutrient_1 : Initial nutrient concentration in the surface layer [mol N/m3]
! Nutrient_2 : Initial nutrient concentration in the bottom layer [mol N/m3]
! bottom_flux : Nutrient replenishment at the bottom [mol N/m3/s]
! nut_relax : Relax nutrient profile to initial profile
! 0: false
! 1: True
! w_s : Plankton settling velocity [m/s] ! w_s : Plankton settling velocity [m/s]
! part_management : algorithm management number particles ! part_management : algorithm management number particles
! mu_r : Respiration rate constant [s-1] ! mu_r : Respiration rate constant [s-1]
...@@ -25,17 +37,23 @@ ...@@ -25,17 +37,23 @@
bio_calc = 1 bio_calc = 1
bio_model = 2 bio_model = 2
mu_max = 2.3148e-5 mu_max = 2.8935e-5
q0 = 5e-02 q0 = 0.05
mortality = 0 mortality = 0
Kn = 9.3870e-6 Kn = 2.14e-4
Vmax = 9.3870e-6 Vmax = 9.3870e-6
Nutrient_0 = 5e-3 nut_profile = 1
nut_layer = 15
Nutrient_0 = 5e-4
Nutrient_1 = 1e-6
Nutrient_2 = 5e-3
bottom_flux = 0
nut_relax = 0
w_s = 0 w_s = 0
part_management = 1 part_management = 0
mu_r = 1.1574e-06 mu_r = 1.1574e-06
kw = 0 kw = 0
q_max = 0.25 q_max = 0.5
m0 = 1e-12 m0 = 1e-12
/ /
...@@ -32,17 +32,17 @@ ...@@ -32,17 +32,17 @@
!------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------
&model_parameter_nml &model_parameter_nml
name_exp ='1depth200' name_exp ='exp3'
depth = 200 depth = 40
dz = 0.5 dz = 0.5
K_value = 1e-6 K_value = 1e-6
time_step = 60 time_step = 60
ndays = 2 ndays = 10
output_time = 2 output_time = 10
mix_state = 1 mix_state = 1
particle_distribution = 2 particle_distribution = 2
nbpart = 5000 nbpart = 20000
z_0 = 10 z_0 = 5
day_length = 16 day_length = 16
rad_max = 600 rad_max = 600
k_bg = 0.4 k_bg = 0.4
......
This diff is collapsed.
...@@ -120,7 +120,7 @@ ...@@ -120,7 +120,7 @@
DN_eul(Nbin)=(K(Nlevel-1)*(N_euler(Nbin-1)-N_euler(Nbin)))/(dz**2) -Nut(Nbin) DN_eul(Nbin)=(K(Nlevel-1)*(N_euler(Nbin-1)-N_euler(Nbin)))/(dz**2) -Nut(Nbin)
N_euler=N_euler+DN_eul*time_step N_euler=N_euler+DN_eul*time_step
N_euler(1:INT(2/dz))=N_euler(1:INT(2/dz))+bottom_flux*time_step
!________________________________________________________________________________ !________________________________________________________________________________
! NUTRIENT CONCENTRATION FOR LAGRANGIAN CALCULATIONS ! NUTRIENT CONCENTRATION FOR LAGRANGIAN CALCULATIONS
...@@ -135,8 +135,19 @@ ...@@ -135,8 +135,19 @@
DN_lag(Nbin)=(K(Nlevel-1)*(N_lag(Nbin-1)-N_lag(Nbin)))/(dz**2) -Nut_taken(Nbin) DN_lag(Nbin)=(K(Nlevel-1)*(N_lag(Nbin-1)-N_lag(Nbin)))/(dz**2) -Nut_taken(Nbin)
N_lag=N_lag+DN_lag*time_step N_lag=N_lag+DN_lag*time_step
N_lag(1:INT(2/dz))=N_lag(1:INT(2/dz))+bottom_flux*time_step
!___________________________________________________________________________________
! RELAXATION
!___________________________________________________________________________________
if (nut_relax.eq.1) then
N_lag=N_init
N_euler=N_init
end if
!================================================================================== !==================================================================================
end subroutine ADV_DIFF end subroutine ADV_DIFF
...@@ -88,6 +88,6 @@ ...@@ -88,6 +88,6 @@
end if end if
end do end do
!=============================================================================== !===============================================================================
end subroutine bio_euler end subroutine bio_euler
...@@ -27,6 +27,7 @@ ...@@ -27,6 +27,7 @@
integer :: i integer :: i
integer :: j integer :: j
integer :: ii
integer :: a integer :: a
integer :: div integer :: div
double precision, allocatable :: dummy1(:) double precision, allocatable :: dummy1(:)
...@@ -36,7 +37,7 @@ ...@@ -36,7 +37,7 @@
double precision :: Pdmax=1 double precision :: Pdmax=1
double precision :: Idark=50 double precision :: Idark=50
double precision, dimension(100000) :: Fi double precision, dimension(100000) :: Fi
double precision :: nb1
allocate(light_level(Nlevel)) allocate(light_level(Nlevel))
allocate(dummy1(Nlevel)) allocate(dummy1(Nlevel))
...@@ -47,7 +48,7 @@ ...@@ -47,7 +48,7 @@
Fi(1:100000)=0 Fi(1:100000)=0
div=0 div=0
nb(1:Nbin)=0 nb(1:Nbin)=0
nb1=0d0
!___________________________________________________________________________ !___________________________________________________________________________
...@@ -85,7 +86,13 @@ ...@@ -85,7 +86,13 @@
weight_cell(i)=weight_cell(i)+Dphyto_lag(i)*time_step*weight_cell(i) !compute growth for each particle weight_cell(i)=weight_cell(i)+Dphyto_lag(i)*time_step*weight_cell(i) !compute growth for each particle
else ! Use Droop's equations else ! Use Droop's equations
do ii=1,20
if ((qlag(i).ge.Qcat(ii)).and.(qlag(i).le.Qcat(ii+1))) then
Q_distrib(bin,ii)=Q_distrib(bin,ii)+1
end if
end do
Uptake_lag(i)=Vmax*(N_lag(bin)/(Kn+N_lag(bin)))*(1-qlag(i)/q_max) ! nutrient uptake Uptake_lag(i)=Vmax*(N_lag(bin)/(Kn+N_lag(bin)))*(1-qlag(i)/q_max) ! nutrient uptake
Nut_taken(bin)=Nut_taken(bin)+(Uptake_lag(i)-kw*qlag(i))*weight(i) ! nutrient taken Nut_taken(bin)=Nut_taken(bin)+(Uptake_lag(i)-kw*qlag(i))*weight(i) ! nutrient taken
if (qlag(i).ge.q0) then !compute growth rate if (qlag(i).ge.q0) then !compute growth rate
...@@ -93,6 +100,12 @@ ...@@ -93,6 +100,12 @@
else else
Dphyto_lag(i)=-mortality -mu_r Dphyto_lag(i)=-mortality -mu_r
end if end if
do ii=1,20
if ((Dphyto_lag(i).ge.growthcat(ii)).and.(Dphyto_lag(i).le.growthcat(ii+1))) then
growth_distrib(bin,ii)=growth_distrib(bin,ii)+1
end if
end do
qlag(i)=qlag(i)+(Uptake_lag(i)-kw*qlag(i)-Dphyto_lag(i)*qlag(i))*time_step !change in cell quota qlag(i)=qlag(i)+(Uptake_lag(i)-kw*qlag(i)-Dphyto_lag(i)*qlag(i))*time_step !change in cell quota
weight_cell(i)=weight_cell(i)+Dphyto_lag(i)*time_step*weight_cell(i) !compute growth for each particle weight_cell(i)=weight_cell(i)+Dphyto_lag(i)*time_step*weight_cell(i) !compute growth for each particle
...@@ -113,7 +126,7 @@ ...@@ -113,7 +126,7 @@
end if ! end if !
end if ! end if !
nb1(bin)=nb1(bin)+1
nb(bin)=nb(bin)+weight(i) nb(bin)=nb(bin)+weight(i)
binned_weight(bin)=binned_weight(bin)+weight(i) binned_weight(bin)=binned_weight(bin)+weight(i)
binned_Qlag(bin)=binned_Qlag(bin) + qlag(i)*weight(i) binned_Qlag(bin)=binned_Qlag(bin) + qlag(i)*weight(i)
...@@ -124,13 +137,15 @@ ...@@ -124,13 +137,15 @@
do i=1,Nbin do i=1,Nbin
Q_distrib(i,:)=Q_distrib(i,:)/nb1(i)
growth_distrib(i,:)=growth_distrib(i,:)/nb1(i)
binned_Qlag(i)=binned_Qlag(i)/nb(i) binned_Qlag(i)=binned_Qlag(i)/nb(i)
binned_growth(i)=binned_growth(i)/nb(i) binned_growth(i)=binned_growth(i)/nb(i)
end do end do
Nut_taken=Nut_taken/dz Nut_taken=Nut_taken/dz
if (part_management.eq.1) then !particle management if (part_management.eq.1) then !particle management
CALL PART_MGMT(div,Fi) CALL PART_MGMT(div,Fi)
......
...@@ -18,8 +18,10 @@ integer :: i,a,b ...@@ -18,8 +18,10 @@ integer :: i,a,b
!======================grid======================================= !======================grid=======================================
time_vec(1)=0d0
do i=2,L_time
time_vec(i)=time_vec(i-1)+(time_step*output_time)/86400d0
end do
levels(1)=0d0 levels(1)=0d0
do i=2,Nlevel do i=2,Nlevel
levels(i)=levels(i-1)+dz levels(i)=levels(i-1)+dz
......
...@@ -42,7 +42,9 @@ ...@@ -42,7 +42,9 @@
if (particle_distribution.eq.1) then ! release all particles in at a particular depth if (particle_distribution.eq.1) then ! release all particles in at a particular depth
z_cell(1:nbpart)=z_0 bin=1+floor(z_0*z_res)
z_cell(1:nbpart)=(levels(bin)+levels(bin+1))/2
else !release all particles uniformely in the domain else !release all particles uniformely in the domain
do i=1,nbpart ! do i=1,nbpart !
z_cell(i)=depth/nbpart/2d0+DBLE(i-1)*(depth/nbpart) ! z_cell(i)=depth/nbpart/2d0+DBLE(i-1)*(depth/nbpart) !
...@@ -61,7 +63,7 @@ ...@@ -61,7 +63,7 @@
weight_cell(1:nbpart)=m0 ! weight_cell(1:nbpart)=m0 !
weight(1:nbpart)=S*weight_cell ! weight(1:nbpart)=S*weight_cell !
binned_weight(1:Nbin)=0
binned_pos=0 binned_pos=0
euler(1:Nbin)=0 euler(1:Nbin)=0
binned_growth(1:Nbin)=0 binned_growth(1:Nbin)=0
...@@ -79,10 +81,30 @@ ...@@ -79,10 +81,30 @@
qlag(1:nbpart)=q0 ! qlag(1:nbpart)=q0 !
binned_Qlag(1:Nbin)=q0 ! binned_Qlag(1:Nbin)=q0 !
q=q*euler ! q=q*euler !
N_euler(1:Nbin)=Nutrient_0 !initial nutrients growthcat(1)=0d0
N_lag(1:Nbin)=Nutrient_0 ! Qcat(1)=q0
Nut_taken(1:Nbin)=0d0 ! do i=2,21
Qcat(i)=Qcat(i-1)+(q_max-q0)/20d0
growthcat(i)=growthcat(i-1)+(mu_max)/20d0
end do
!_________________________________________________________________________________
!nutrient profile
!_________________________________________________________________________________
if (nut_profile.eq.1) then
N_init(1:Nbin)=Nutrient_0
N_init(1:Nbin)=Nutrient_0
else
N_init(1:INT((depth-nut_layer)/dz))=Nutrient_2
N_init(INT((depth-nut_layer)/dz):Nbin)=Nutrient_1
end if
N_lag=N_init
N_euler=N_init
Nut_taken(1:Nbin)=0d0
!============================================================================== !==============================================================================
......
...@@ -37,6 +37,8 @@ ...@@ -37,6 +37,8 @@
Qeuler=trim(root)//'/EULER_CELL_QUOTA.dat' !eulerian cell quota Qeuler=trim(root)//'/EULER_CELL_QUOTA.dat' !eulerian cell quota
Qlagrange=trim(root)//'/LAGRANGE_CELL_QUOTA.dat' !lagrangian cell quota Qlagrange=trim(root)//'/LAGRANGE_CELL_QUOTA.dat' !lagrangian cell quota
turb=trim(root)//'/DIFFUSIVITY.dat' !diffusivity profile turb=trim(root)//'/DIFFUSIVITY.dat' !diffusivity profile
depth_vector=trim(root)//'/DEPTH.dat' !time vector
time_vector=trim(root)//'/TIME.dat' !depth vector
open(12,file=conc_euler) !open output files open(12,file=conc_euler) !open output files
...@@ -49,6 +51,8 @@ ...@@ -49,6 +51,8 @@
open(32,file=Qlagrange) ! open(32,file=Qlagrange) !
open(33,file=growth_eul) ! open(33,file=growth_eul) !
open(34,file=growth_lag) ! open(34,file=growth_lag) !
open(35,file=time_vector)
open(36,file=depth_vector)
!______________________________________________________________________________ !______________________________________________________________________________
...@@ -99,8 +103,8 @@ ...@@ -99,8 +103,8 @@
write(32,*) binned_Qlag ! write(32,*) binned_Qlag !
write(33,*) grow ! write(33,*) grow !
write(34,*) binned_growth ! write(34,*) binned_growth !
Q_distrib2(:,:,1)=0d0
i=2
do time=1,nsteps !START THE TIME LOOP! do time=1,nsteps !START THE TIME LOOP!
...@@ -109,13 +113,14 @@ ...@@ -109,13 +113,14 @@
end if ! end if !
CALL daylight !compute the incident light CALL daylight !compute the incident light
CALL rdm_walk !perform the random walk CALL rdm_walk !perform the random walk
CALL ADV_DIFF !compute advection and diffusion in a eulerian way
if (bio_calc.eq.1) then ! if (bio_calc.eq.1) then !
CALL light_growth !lagrangian growth CALL light_growth !lagrangian growth
CALL bio_euler !eulerian growth CALL bio_euler !eulerian growth
end if end if
CALL ADV_DIFF !compute advection and diffusion in a eulerian way
if (mod(time,output_time).eq.0) then !write outputs each output_time if (mod(time,output_time).eq.0) then !write outputs each output_time
! !
write(21,*) K ! write(21,*) K !
...@@ -129,6 +134,9 @@ ...@@ -129,6 +134,9 @@
if(bio_model.eq.2) then ! if(bio_model.eq.2) then !
write(32,*) binned_Qlag ! write(32,*) binned_Qlag !
write(31,*) q/euler ! write(31,*) q/euler !
Q_distrib2(:,:,i)=Q_distrib
growth_distrib2(:,:,i)=growth_distrib
i=i+1
end if ! end if !
! !
end if ! end if !
...@@ -140,14 +148,17 @@ ...@@ -140,14 +148,17 @@
binned_Qlag(1:Nbin)=0d0 ! binned_Qlag(1:Nbin)=0d0 !
binned_growth(1:Nbin)=0d0 ! binned_growth(1:Nbin)=0d0 !
binned_weight(1:Nbin)=0d0 ! binned_weight(1:Nbin)=0d0 !
Q_distrib=0d0
growth_distrib=0d0
CALL clock !display the time of the simulation CALL clock !display the time of the simulation
end do !END OF THE TIME LOOP! end do !END OF THE TIME LOOP!
write(35,*)time_vec
write(36,*)levels(1:Nlevel-1)
CALL write_output
close(21) ! close output files close(21) ! close output files
close(13) ! close(13) !
......
...@@ -31,6 +31,8 @@ ...@@ -31,6 +31,8 @@
integer :: day !for the clock integer :: day !for the clock
integer :: output_time=100 !number of step for output integer :: output_time=100 !number of step for output
integer :: minute !for the clock integer :: minute !for the clock
integer :: L_time !length of output file
double precision, allocatable :: time_vec(:) !time vector
!_________________________________________________________________________________ !_________________________________________________________________________________
! SPACE PARAMETERS ! SPACE PARAMETERS
...@@ -63,6 +65,9 @@ ...@@ -63,6 +65,9 @@
character(len=400) :: Qlagrange character(len=400) :: Qlagrange
character(len=400) :: growth_eul character(len=400) :: growth_eul
character(len=400) :: growth_lag character(len=400) :: growth_lag
character(len=400) :: time_vector
character(len=400) :: depth_vector
!________________________________________________________________________________ !________________________________________________________________________________
...@@ -107,6 +112,13 @@ ...@@ -107,6 +112,13 @@
double precision :: w_s=2.3148e-06 double precision :: w_s=2.3148e-06
double precision :: mortality=0.00006944 double precision :: mortality=0.00006944
double precision :: Nutrient_0=6 double precision :: Nutrient_0=6
double precision :: Nutrient_1=2
double precision :: Nutrient_2=2
integer :: nut_profile=1
integer :: nut_relax=1
double precision :: nut_layer=10
double precision :: bottom_flux=1e-6
double precision, allocatable :: N_init(:)
double precision, allocatable :: weight(:) double precision, allocatable :: weight(:)
double precision, allocatable :: binned_weight(:) double precision, allocatable :: binned_weight(:)
double precision, allocatable :: Uptake_eul(:) double precision, allocatable :: Uptake_eul(:)
...@@ -137,6 +149,12 @@ ...@@ -137,6 +149,12 @@
double precision :: q_max=5e-03 double precision :: q_max=5e-03
double precision, allocatable :: S(:) double precision, allocatable :: S(:)
double precision, allocatable :: weight_cell(:) double precision, allocatable :: weight_cell(:)
double precision, allocatable :: Q_distrib(:,:)
double precision, allocatable :: Q_distrib2(:,:,:)
double precision, dimension(21) :: Qcat
double precision, dimension(21) :: growthcat
double precision, allocatable :: growth_distrib(:,:)
double precision, allocatable :: growth_distrib2(:,:,:)
contains contains
...@@ -152,8 +170,9 @@ ...@@ -152,8 +170,9 @@
nbpart,z_0,day_length,rad_max,k_bg nbpart,z_0,day_length,rad_max,k_bg
namelist /bio_parameter_nml/ bio_calc,bio_model,mu_max,q0,mortality, & namelist /bio_parameter_nml/ bio_calc,bio_model,mu_max,q0,mortality, &
Kn,Vmax,Nutrient_0,w_s,part_management,mu_r, & Kn,Vmax,nut_profile,nut_layer,Nutrient_0, &
kw,q_max,m0 Nutrient_1,Nutrient_2,bottom_flux,nut_relax, &
w_s,part_management,mu_r,kw,q_max,m0
open(14,file='nml/model_parameter.nml',status='old') open(14,file='nml/model_parameter.nml',status='old')
read(14,nml=model_parameter_nml) read(14,nml=model_parameter_nml)
...@@ -184,8 +203,8 @@ ...@@ -184,8 +203,8 @@
Nlevel=INT(depth/dz)+1 Nlevel=INT(depth/dz)+1
Nbin=Nlevel-1 Nbin=Nlevel-1
z_res=INT(Nbin/depth) z_res=INT(Nbin/depth)
L_time=INT(nsteps/output_time)+1
allocate(levels(Nlevel)) allocate(levels(Nlevel))
allocate(weight(nbpart)) allocate(weight(nbpart))
allocate(z_cell(nbpart)) allocate(z_cell(nbpart))
...@@ -196,6 +215,7 @@ ...@@ -196,6 +215,7 @@
allocate(K(Nlevel)) allocate(K(Nlevel))
allocate(Uptake_eul(Nbin)) allocate(Uptake_eul(Nbin))
allocate(q(Nbin)) allocate(q(Nbin))
allocate(N_init(Nbin))
allocate(N_euler(Nbin)) allocate(N_euler(Nbin))
allocate(N_lag(Nbin)) allocate(N_lag(Nbin))
allocate(Uptake_lag(nbpart)) allocate(Uptake_lag(nbpart))
...@@ -210,7 +230,11 @@ ...@@ -210,7 +230,11 @@
allocate(grow(Nbin)) allocate(grow(Nbin))
allocate(Nut(Nbin)) allocate(Nut(Nbin))
allocate(binned_growth(Nbin)) allocate(binned_growth(Nbin))
allocate(time_vec(L_time))
allocate(Q_distrib(Nbin,20))
allocate(Q_distrib2(Nbin,20,L_time))
allocate(growth_distrib(Nbin,20))
allocate(growth_distrib2(Nbin,20,L_time))
end subroutine end subroutine
......
!________________________________________________________________________________
!DESCRIPTION: This routine create the NETCDF output
!file. You must have the FORTRAN NETCDF libraries
!installed on your computer to use this routine.
!________________________________________________________________________________
!INTERFACE:
subroutine write_output
!MODULE USES:
use netcdf