Skip to content

GitLab

Commit b6146951 authored by root's avatar root
Browse files

nouveau commit

parent 8a16b803
Hide whitespace changes
Inline Side-by-side
nml/bio_parameter.nml
View file @ b6146951
......@@ -36,24 +36,25 @@
&bio_parameter_nml
bio_calc = 1
bio_model = 2
mu_max = 2.8935e-5
photo_model = 2
bio_model = 1
mu_max = 2.3148e-5
q0 = 0.05
mortality = 0
Kn = 2.14e-4
Vmax = 9.3870e-6
Vmax = 1.1574e-5
nut_profile = 1
nut_layer = 15
Nutrient_0 = 5e-4
Nutrient_0 = 5e-3
Nutrient_1 = 1e-6
Nutrient_2 = 5e-3
bottom_flux = 0
nut_relax = 0
w_s = 0
w_s = 2.3148e-06
part_management = 0
mu_r = 1.1574e-06
kw = 0
q_max = 0.5
m0 = 1e-12
mu_r = 1.1574e-6
kw = 1.1574e-7
q_max = 0.25
m0 = 4e-6
/
nml/model_parameter.nml
View file @ b6146951
......@@ -32,12 +32,12 @@
!------------------------------------------------------------------------------------------
&model_parameter_nml
name_exp ='exp3'
name_exp ='test_acclim'
depth = 40
dz = 0.5
K_value = 1e-6
time_step = 60
ndays = 10
ndays = 15
output_time = 10
mix_state = 1
particle_distribution = 2
......
src/ADV_DIFF.f90
View file @ b6146951
......@@ -27,21 +27,49 @@
double precision, allocatable :: DC(:)
double precision, allocatable :: DN_eul(:)
double precision, allocatable :: DN_lag(:)
double precision, allocatable :: Dq(:)
double precision, allocatable :: q2(:)
integer :: i
double precision, allocatable :: Dq(:),DQQ(:,:)
double precision, allocatable :: q2(:),DNcell(:)
integer :: i,j
double precision :: alpha
double precision :: beta
double precision :: teta
double precision, allocatable :: DCL(:)
Allocate(DC(Nbin))
Allocate(DN_eul(Nbin))
Allocate(DN_lag(Nbin))
Allocate(Dq(Nbin))
Allocate(q2(Nbin))
Allocate(DQQ(Nbin,Ncat))
allocate(DNcell(Nbin))
Allocate(DCL(Nbin))
if(photo_model.eq.2) then
do i=2,Nbin-1
DCL(i)=w_s*(CL(i+1)-CL(i))/dz
end do
DCL(1)=w_s*CL(2)/dz
DCL(Nbin)=-w_s*CL(Nbin)/dz
CL=CL+DCL*time_step
do i=2,Nbin-1
DCL(i)=(K(i)*(CL(i-1)-CL(i))+K(i+1)*(CL(i+1)-CL(i)))/(dz**2)
end do
DCL(1)=(K(2)*(CL(2)-CL(1)))/(dz**2)
DCL(Nbin)=(K(Nlevel-1)*(CL(Nbin-1)-CL(Nbin)))/(dz**2)
CL=CL+DCL*time_step
end if
if(bio_model.le.2) then
!________________________________________________________________________________
! ADVECTION OF CELL QUOTA FOR EULERIAN CALCULATIONS
!________________________________________________________________________________
......@@ -49,13 +77,13 @@
do i=2,Nbin-1
Dq(i)=w_s*(q(i+1)-q(i))/dz
DNcell(i)=w_s*(Ncell(i+1)-Ncell(i))/dz
end do
Dq(1)=w_s*q(2)/dz
Dq(Nbin)=-w_s*q(Nbin)/dz
DNcell(1)=w_s*Ncell(2)/dz
DNcell(Nbin)=-w_s*Ncell(Nbin)/dz
q=q+Dq*time_step
Ncell=Ncell+DNcell*time_step
!________________________________________________________________________________
......@@ -65,13 +93,13 @@
do i=2,Nbin-1
Dq(i)=(K(i)*(q(i-1)-q(i))+K(i+1)*(q(i+1)-q(i)))/(dz**2)
DNcell(i)=(K(i)*(Ncell(i-1)-Ncell(i))+K(i+1)*(Ncell(i+1)-Ncell(i)))/(dz**2)
end do
Dq(1)=(K(2)*(q(2)-q(1)))/(dz**2)
Dq(Nbin)=(K(Nlevel-1)*(q(Nbin-1)-q(Nbin)))/(dz**2)
DNcell(1)=(K(2)*(Ncell(2)-Ncell(1)))/(dz**2)
DNcell(Nbin)=(K(Nlevel-1)*(Ncell(Nbin-1)-Ncell(Nbin)))/(dz**2)
q=q+Dq*time_step
Ncell=Ncell+DNcell*time_step
......@@ -106,18 +134,33 @@
euler=euler+DC*time_step
else
do j=1,Ncat
do i=2,Nbin-1
DQQ(i,j)=(K(i)*(q_dist(i-1,j)-q_dist(i,j))+K(i+1)*(q_dist(i+1,j)-q_dist(i,j)))/(dz**2)
end do
DQQ(1,j)=(K(2)*(q_dist(2,j)-q_dist(1,j)))/(dz**2)
DQQ(Nbin,j)=(K(Nlevel-1)*(q_dist(Nbin-1,j)-q_dist(Nbin,j)))/(dz**2)
q_dist(:,j)=q_dist(:,j)+DQQ(:,j)*time_step
end do
end if
!________________________________________________________________________________
! NUTRIENT CONCENTRATION FOR EULERIAN CALCULATIONS
!________________________________________________________________________________
do i=2,Nbin-1
DN_eul(i)=(K(i)*(N_euler(i-1)-N_euler(i))+K(i+1)*(N_euler(i+1)-N_euler(i)))/(dz**2) -Nut(i)
DN_eul(i)=(K(i)*(N_euler(i-1)-N_euler(i))+K(i+1)*(N_euler(i+1)-N_euler(i)))/(dz**2)-Nut(i)
end do
DN_eul(1)=(K(2)*(N_euler(2)-N_euler(1)))/(dz**2) -Nut(1)
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(1:INT(2/dz))=N_euler(1:INT(2/dz))+bottom_flux*time_step
......
src/BIO_EULER.f90
View file @ b6146951
......@@ -22,19 +22,30 @@
IMPLICIT NONE
integer :: i
integer :: i,j,jj
integer :: a
double precision, allocatable :: light_level(:)
double precision, allocatable :: light_bin(:)
double precision, allocatable :: dummy1(:)
double precision, allocatable :: dummy1(:),q_dummy(:),mu_dummy(:),nt(:),dt1(:)
double precision :: euler_av
double precision :: Pdmax=1
double precision :: Idark=50
double precision, allocatable :: Ndummy(:),Pl(:),Pd(:),CLQUOT(:)
allocate(light_level(Nlevel))
allocate(dummy1(Nlevel))
allocate(light_bin(Nbin))
allocate(q_dummy(Ncat))
allocate(mu_dummy(Ncat))
allocate(nt(Ncat))
allocate(dt1(Ncat))
allocate(Ndummy(Nbin))
allocate(Pl(Nbin))
allocate(Pd(Nbin))
allocate(CLQUOT(Nbin))
!____________________________________________________________________________
......@@ -55,39 +66,142 @@
do i=1,Nbin !compute light intensity in each grid bin
light_bin(i)=(light_level(i)+light_level(i+1))/2 !
end do !
!_____________________________________________________________________________
! GROWTH
!_____________________________________________________________________________
Nut=0d0
mu_dist=0d0
do i=1,Nbin
if (photo_model.eq.1) then
growlight(i)=Pdmax*(1d0-dexp(-light_bin(i)/Idark))
else
CLQUOT=CL/euler
!growlight(i)=Pdmax*(1d0-dexp(-light_bin(i)/(Idark+1000*CLQUOT(i))))
growlight(i)=light_bin(i)/(light_bin(i)+(500*CLQUOT(i))**2)
if (light_bin(i)/(rad_max*0.45*4.16).lt.0.01) then
Clquot_ACC(i)=0.04
else
Clquot_ACC(i)=0.01+0.03*log(light_bin(i)/(rad_max*0.45*4.16))/log(0.01)
end if
growlight(i)=Pdmax*(1d0-dexp(-light_bin(i)/Idark)) !calculate P(I)
end if
q(i)=Ncell(i)/euler(i)
if(q(i).gt.q_max) then
!write(*,*)q(i)
Ncell(i)=q_max*euler(i)
q(i)=Ncell(i)/euler(i)
end if
if (bio_model.eq.1) then ! Monod equations
grow(i)=(growlight(i)*(mu_max*(N_euler(i))/(Kn+N_euler(i))) -mortality) -mu_r ! growth rate
Nut(i)=(growlight(i)*mu_max*N_euler(i)/(Kn+N_euler(i)))*q0*euler(i) ! nutrients taken
Nut(i)=(growlight(i)*mu_max*(N_euler(i)/(Kn+N_euler(i))))*q0*euler(i)
euler(i)=euler(i)+grow(i)*time_step*euler(i) !compute the growth with euler solver
else !Droop equations
CL(i)=CL(i)+(1/tau*(Clquot_ACC(i)-CLQUOT(i))*euler(i)+grow(i)*CL(i))*time_step
elseif (bio_model.eq.4) then
CLQUOT_dumm=0
do j=1,nCLcat
if ((Clquot_ACC(i).gt.CLQUOTDIST(i,j).and.(j.ne.nCLcat).and.(CLQUOTDIST(i,j).gt.0d0)))then
CLQUOT_dumm(j+1)=CLQUOT_dumm(j+1)+(1/tau*(Clquot_ACC(i)-CLQUOTDIST(i,j))/(CLquat(i,j+1)-CLquat(i,j)))*time_step
elseif ((Clquot_ACC(i).lt.CLQUOTDIST(i,j).and.(j.ne.1).and.(CLQUOTDIST(i,j).gt.0d0)) then
CLQUOT_dumm(j-1)=CLQUOT_dumm(j-1)+(1/tau*(Clquot_ACC(i)-CLQUOTDIST(i,j))/(CLquat(i,j)-CLquat(i,j-1)))*time_step
end if
end do
CL_to_mu=(light_bin(i)/(light_bin(i)+(500*CLcat)**2))*(mu_max*(N_euler(i))/(Kn+N_euler(i)))-mortality -mu_r
Nut(i)=sum(((light_bin(i)/(light_bin(i)+(500*CLcat)**2))*mu_max*(N_euler(i)/(Kn+N_euler(i))))*q0*CLQUOTDIST(i,:))
CLQUOTDIST(i,:)=CLQUOTDIST(i,:)+CLQUOT_dumm
CLQUOTDIST(i,:)=CLQUOTDIST(i,:)+CL_to_mu*CLQUOTDIST(i,:)*time_step
Uptake_eul(i)=Vmax*N_euler(i)/(Kn+N_euler(i))*(1-(q(i)/euler(i))/q_max) !nutrient uptake
q(i)=q(i)+(Uptake_eul(i)-kw*q(i))*time_step*euler(i) !change in cell quota
Nut(i)=(Uptake_eul(i)-kw*q(i))*euler(i) !nutrient taken
euler(i)=sum(CLQUOTDIST(i,:))
if (q(i).ge.q0*euler(i)) then !compute the growth rate
grow(i)=growlight(i)*mu_max*(1-q0/(q(i)/euler(i)))-mu_r !
elseif (bio_model.eq.2) then !Droop equations
Uptake_eul(i)=Vmax*N_euler(i)/(Kn+N_euler(i))*(1-q(i)/q_max)
Ncell(i)=Ncell(i)+(Uptake_eul(i)-kw*q(i))*time_step*euler(i) -mu_r*Ncell(i)*time_step
Nut(i)=(Uptake_eul(i)-kw*q(i))*euler(i)-mu_r*Ncell(i) !nutrient taken
if (q(i).gt.q0) then !compute the growth rate
grow(i)=growlight(i)*mu_max*(1-q0/q(i))-mu_r
!grow(i)=growlight(i)*mu_max*(q(i)/q_max)-mu_r
else !
grow(i)=-mortality -mu_r !
end if !
end if !
euler(i)=euler(i)+grow(i)*time_step*euler(i) !compute the growth with euler solver
elseif (bio_model.eq.3) then
q_dummy=0
nt=0
dt1=0
do j=1,Ncat
uptake_dist(j)=Vmax*N_euler(i)/(Kn+N_euler(i))*(1-q_cat(j)/q_max)
Nut(i)=Nut(i)+uptake_dist(j)*q_dist(i,j)-kw*q_cat(j)*q_dist(i,j)
if (q_cat(j).gt.q0) then
q_to_mu(j)=growlight(i)*mu_max*(1-q0/q_cat(j))-mu_r
else
q_to_mu(j)=-mortality -mu_r
end if
q_to_mu(1)=-mortality -mu_r
if((uptake_dist(j).lt.q_to_mu(j)*q_cat(j)+kw*q_cat(j)).and.(j.ne.1).and.(q_dist(i,j).gt.0d0)) then
dt1(j)=(q_cat(j-1)-q_cat(j))/(uptake_dist(j)-q_to_mu(j)*q_cat(j)-kw*q_cat(j))
nt(j)=dt1(j)/time_step;
q_dummy(j-1)=q_dummy(j-1)+q_dist(i,j)/nt(j)
q_dummy(j)=q_dummy(j)-q_dist(i,j)/nt(j)
elseif((uptake_dist(j).gt.q_to_mu(j)*q_cat(j)+kw*q_cat(j)).and.(j.ne.Ncat).and.(q_dist(i,j).gt.0d0)) then
dt1(j)=(q_cat(j+1)-q_cat(j))/(uptake_dist(j)-q_to_mu(j)*q_cat(j)-kw*q_cat(j))
nt(j)=dt1(j)/time_step;
q_dummy(j+1)=q_dummy(j+1)+q_dist(i,j)/nt(j)
q_dummy(j)=q_dummy(j)-q_dist(i,j)/nt(j)
end if
end do
! do j=1,Ncat
! do jj=1,Ncat
! if((q_to_mu(j).ge.mu_edge(jj)).and.(q_to_mu(j).le.mu_edge(jj+1))) then
! mu_dist(i,jj)=mu_dist(i,jj)+q_dist(i,j)
! end if
! end do
! end do
! write(*,*)euler(i),sum(q_dist(i,:))+ sum(q_dist(i,:)*q_to_mu*time_step)
q_dist(i,:)=q_dist(i,:)+ q_dist(i,:)*q_to_mu*time_step
q_dist(i,:)=q_dist(i,:)+q_dummy
euler(i)=sum(q_dist(i,:))
end if
end do
!===============================================================================
end subroutine bio_euler
src/BIO_LAGRANGE.f90
View file @ b6146951
......@@ -38,12 +38,18 @@
double precision :: Idark=50
double precision, dimension(100000) :: Fi
double precision,allocatable :: nb1(:)
double precision,allocatable :: Pl(:),Pd(:)
double precision, allocatable :: CLQUOT_lag(:)
allocate(light_level(Nlevel))
allocate(dummy1(Nlevel))
allocate(light_bin(Nbin))
allocate(nb(Nbin))
allocate(nb1(Nbin))
allocate(Pl(Nbpart))
allocate(Pd(Nbpart))
allocate(CLQUOT_lag(Nbpart))
Fi(1:100000)=0
div=0
......@@ -76,51 +82,87 @@
do i=1,nbpart
bin=1+floor(z_cell(i)*z_res) !find the bin where the particle is.
qlag(i)=Ncell_lag(i)/weight(i)
if (qlag(i).gt.q_max)then
!write(*,*)qlag(i)
Ncell_lag(i)=q_max*weight(i)
qlag(i)=Ncell_lag(i)/weight(i)
end if
bin=1+floor(z_cell(i)*z_res)
!find the bin where the particle is.
if (photo_model.eq.1) then
uptlight(i)=Pdmax*(1d0-dexp(-light_bin(bin)/Idark)) !compute P(I)
else
CLQUOT_lag(i)=CL_lag(i)/weight(i)
! uptlight(i)=Pdmax*(1d0-dexp(-light_bin(bin)/(Idark+1000*CLQUOT_lag(i))))
uptlight(i)=light_bin(bin)/(light_bin(bin)+(500*CLQUOT_lag(i))**2)
if (light_bin(bin)/(rad_max*0.45*4.16).lt.0.01) then
Clquot_ACC(bin)=0.04
else
Clquot_ACC(bin)=0.01+0.03*log(light_bin(bin)/(rad_max*0.45*4.16))/log(0.01)
end if
end if
if (bio_model.eq.1) then ! Use Monod's equations
Dphyto_lag(i)=(uptlight(i)*mu_max*(N_lag(bin)/(Kn+N_lag(bin))) -mortality)-mu_r !compute growth rate
Nut_taken(bin)=Nut_taken(bin)+(uptlight(i)*mu_max*(N_lag(bin)/(Kn+N_lag(bin))))*q0*weight(i) !nutrient taken
weight_cell(i)=weight_cell(i)+Dphyto_lag(i)*time_step*weight_cell(i) !compute growth for each particle
weight(i)=weight(i)+Dphyto_lag(i)*time_step*weight(i) !compute growth for each particle
if (photo_model.eq.2) then
CL_lag(i)=CL_lag(i)+(1/tau*(Clquot_ACC(bin)-CLQUOT_lag(i))*weight(i)+Dphyto_lag(i)*CL_lag(i))*time_step
end if
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
! do ii=1,Ncat
! if ((qlag(i).ge.q_edge(ii)).and.(qlag(i).le.q_edge(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
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
Ncell_lag(i)=Ncell_lag(i)+(Uptake_lag(i)-kw*qlag(i))*time_step*weight(i)-mu_r*Ncell_lag(i)*time_step
Nut_taken(bin)=Nut_taken(bin)+(Uptake_lag(i)-kw*qlag(i))*weight(i)-mu_r*Ncell_lag(i) ! nutrient taken
if (qlag(i).ge.q0) then
Dphyto_lag(i)=uptlight(i)*mu_max*(1-q0/qlag(i)) -mu_r
!Dphyto_lag(i)=uptlight(i)*mu_max*(qlag(i)/q_max) -mu_r
else
Dphyto_lag(i)=-mortality -mu_r
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
weight_cell(i)=weight_cell(i)+Dphyto_lag(i)*time_step*weight_cell(i) !compute growth for each particle
! do ii=1,Ncat
!if ((Dphyto_lag(i).ge.mu_edge(ii)).and.(Dphyto_lag(i).le.mu_edge(ii+1))) then
! growth_distrib(bin,ii)=growth_distrib(bin,ii)+1
! end if
! end do
weight(i)=weight(i)+Dphyto_lag(i)*time_step*weight(i) !compute growth for each particle
end if
if(weight_cell(i).gt.2d0*m0) then !cell division
S(i)=2*S(i) !
weight_cell(i)=weight_cell(i)*0.5 !
end if !
weight(i)=S(i)*weight_cell(i) !
if (part_management.eq.1) then !particle management
if(S(i).gt.16*4e6) then !if 4 divisions
if(weight(i).gt.16*m0) then !if 4 divisions
div=div+1 ! Number of particle to be divided
Fi(div)=i ! save which particle to be divided
end if !
......
src/INITIALIZATION.f90
View file @ b6146951
......@@ -23,7 +23,14 @@
integer :: i
integer :: j
integer :: a
double precision, allocatable ::light_level(:),light_bin(:)
double precision, allocatable :: dummy1(:)
allocate(light_level(Nlevel))
allocate(light_bin(Nbin))
allocate(dummy1(Nlevel))
!_____________________________________________________________________________
......@@ -50,7 +57,7 @@
z_cell(i)=depth/nbpart/2d0+DBLE(i-1)*(depth/nbpart) !
end do !
end if !
!______________________________________________________________________________
......@@ -59,9 +66,11 @@
S(1:nbpart)=4e6 !initial biomass per particle
weight_cell(1:nbpart)=m0 !
weight(1:nbpart)=S*weight_cell !
!S(1:nbpart)=4e6 !initial biomass per particle
!weight_cell(1:nbpart)=m0 !
!weight(1:nbpart)=S*weight_cell !
weight(1:nbpart)=m0
binned_weight(1:Nbin)=0
binned_pos=0
......@@ -77,18 +86,55 @@
euler=binned_weight/dz
q(1:Nbin)=q0 !initial cell quota
qlag(1:nbpart)=q0 !
Ncell(1:Nbin)=q0*euler !initial cell quota
Ncell_lag=q0*weight
qlag(1:nbpart)=Ncell_lag/weight !
binned_Qlag(1:Nbin)=q0 !
q=q*euler !
q=Ncell/euler !
growthcat(1)=0d0
Qcat(1)=q0
do i=2,21
Qcat(i)=Qcat(i-1)+(q_max-q0)/20d0
growthcat(i)=growthcat(i-1)+(mu_max)/20d0
mu_edge(1)=-2e-6
q_edge(1)=q0
do i=2,Ncat+1
q_edge(i)=q_edge(i-1)+(q_max-q0)/Ncat
mu_edge(i)=mu_edge(i-1)+(mu_max+2e-6)/Ncat
end do
do i=1,Ncat
q_cat(i)=(q_edge(i)+q_edge(i+1))/2d0
mu_cat(i)=(mu_edge(i)+mu_edge(i+1))/2d0
end do
q_dist=0d0
q_dist(:,1)=euler(10)
mu_dist=0
if (photo_model.eq.2) then
light_level=rad*0.45*4.16*dexp(-k_bg*levels)
a=0
do i=1,Nlevel
dummy1(Nlevel-a)=light_level(i)
a=a+1
end do !
light_level=dummy1
do i=1,Nbin
light_bin(i)=(light_level(i)+light_level(i+1))/2
if (light_bin(i)/(rad_max*0.45*4.16).lt.0.01) then