Skip to content

Commit

Permalink
physiologycal processes improved
Browse files Browse the repository at this point in the history
  • Loading branch information
@zhenwu0728
zhenwu0728 committed Jun 1, 2019
1 parent 5ef6259 commit dd045d3
Show file tree
Hide file tree
Showing 7 changed files with 212 additions and 117 deletions.
2 changes: 1 addition & 1 deletion agent_div.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,7 @@ function agent_move(phyts_a,velᵈ,g,ΔT::Int64)
# phyt.y = max(1.5,min(g.Ny-0.5,phyt.y - dy*(1+rand()/5)))
phyt.x = phyt.x - dx*(1+rand()/3)
phyt.y = phyt.y - dy*(1+rand()/3)
phyt.z = max(1.0,min(g.Nz-0.1,phyt.z - dz*(1+rand()/3)))
phyt.z = max(2.0,min(g.Nz-0.1,phyt.z - dz*(1+rand()/3)))
# periodic domian
if phyt.x ≥ g.Nx - 0.5
phyt.x = phyt.x - g.Nx + 2.0
Expand Down
24 changes: 12 additions & 12 deletions model_setup.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,33 +3,33 @@ function setup_agents(N::Int64,Cquota::Array,Nn::Int64,mean::Float64,var::Float6
phyts0 = DataFrame(x=Float64[], y=Float64[], z=Float64[], gen=Int64[], size=Float64[], Cq1=Float64[], Cq2=Float64[], Nq=Float64[], chl=Float64[], sp=Int64[])
for i in 1:N
# agent location
x = rand(40*grid.Nx:60*grid.Nx)/100
y = rand(40*grid.Nx:grid.Ny*50)/100
x = rand(30*grid.Nx:70*grid.Nx)/100
y = rand(30*grid.Nx:grid.Ny*50)/100
z = rand(3.5*10:grid.Nz*8)/10
# a normal distribution with mean variance
radm = max(0.05, rand(Normal(mean,var)))
gen = 1
size = radm
Cq1 = Cquota[1]*Nn # Nn is the number of cells one super agent repersents
Cq2 = Cquota[1]*Nn*radm
Nq = 13/106*2*Cq2*Nn
chl = Cq2*0.4*Nn
Nq = 13/106*2*Cq2
chl = Cq2*0.4
sp = 1
push!(phyts0,(x=x,y=y,z=z,gen=gen,size=size,Cq1=Cq1,Cq2=Cq2,Nq=Nq,chl=chl,sp=sp))
end
for i in N+1:2N
# agent location
x = rand(40*grid.Nx:60*grid.Nx)/100
y = rand(50*grid.Ny:60*grid.Ny)/100
x = rand(30*grid.Nx:70*grid.Nx)/100
y = rand(50*grid.Ny:70*grid.Ny)/100
z = rand(3.5*10:grid.Nz*8)/10
# a normal distribution with mean variance
radm = max(0.05, rand(Normal(mean,var)))
gen = 1
size = radm
Cq1 = Cquota[2]*Nn
Cq2 = Cquota[2]*Nn*radm
Nq = 13/106*2*Cq2*Nn
chl = Cq2*0.4*Nn
Nq = 13/106*2*Cq2
chl = Cq2*0.4
sp = 2
push!(phyts0,(x=x,y=y,z=z,gen=gen,size=size,Cq1=Cq1,Cq2=Cq2,Nq=Nq,chl=chl,sp=sp))
end
Expand All @@ -46,17 +46,17 @@ function setup_nutrients(g,nut)
DON = zeros(g.Nx, g.Ny, g.Nz)
POC = zeros(g.Nx, g.Ny, g.Nz)
PON = zeros(g.Nx, g.Ny, g.Nz)
for i in trunc(Int,0.3*g.Nx):trunc(Int,0.7*g.Nx)
for j in trunc(Int,0.3*g.Ny):trunc(Int,0.7*g.Ny)
for k in 1:Int(0.6*g.Nz)
for i in trunc(Int,0.25*g.Nx):trunc(Int,0.75*g.Nx)
for j in trunc(Int,0.25*g.Ny):trunc(Int,0.75*g.Ny)
for k in 1:Int(0.625*g.Nz)
DIC[i, j, k] = DIC[i, j, k] + nut[1]
DIN[i, j, k] = DIN[i, j, k] + nut[2]*0.5
DOC[i, j, k] = DOC[i, j, k] + nut[3]
DON[i, j, k] = DON[i, j, k] + nut[4]
POC[i, j, k] = POC[i, j, k] + nut[5]
PON[i, j, k] = PON[i, j, k] + nut[6]
end
for k in Int(0.6*g.Nz)+1:Int(0.80*g.Nz)
for k in Int(0.625*g.Nz)+1:Int(0.80*g.Nz)
DIC[i, j, k] = DIC[i, j, k] + nut[1]
DIN[i, j, k] = DIN[i, j, k] + nut[2]*1.0
DOC[i, j, k] = DOC[i, j, k] + nut[3]
Expand Down
99 changes: 91 additions & 8 deletions model_update.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,30 +58,31 @@ itList = collect(itvalLo:144:itvalHi);
tN = 4056; # starting time
vfroot = "/nobackup1b/users/jahn/hinpac/grazsame3/run/run.0354/offline-0604/"; # directory of velocity fields

N = 80000 # Number of initial individuals of each species
N = 100000 # Number of initial individuals of each species
Nsp = 2 # Number of species
Nn = Int(1e10) # Number of cells one super-agent represents
Nn = Int(1e15) # Number of cells one super-agent represents
B=setup_agents(N,Cquota,Nn,1.1,0.18,g) # Normal distribution with mean and variance
# model initialization
# create output file
output = create_output(B);
nut = [2.0, 0.1, 20.0, 0.0, 0.0, 0.0] #DIC, DIN, DOC, DON, POC, PON, mmol/m3
nut = [2.0, 0.05, 20.0, 0.0, 0.0, 0.0] #DIC, DIN, DOC, DON, POC, PON, mmol/m3
nutrients= setup_nutrients(g,nut)
remin = rem(kDOC,kDON,kPOC,kPON)
remin = rem(kDOC,kDON,kPOC,kPON);

for t in 1:nTime
phyts_a = copy(B[t]) # read data from last time step
phyts_b,dvid_ct,graz_ct,death_ct,consume=phyt_update(t, ΔT, g, phyts_a, nutrients, IR, temp)
velᵇ = read_offline_vels(vfroot,itList,tN,trunc(Int,t*ΔT/3600));
velᵈ = double_grid(velᵇ,g)
agent_move(phyts_a,velᵈ,g,ΔT)
phyts_b,dvid_ct,graz_ct,consume=phyt_update(t, ΔT, g, phyts_a, nutrients, IR, temp)
agent_move(phyts_b,velᵈ,g,ΔT)
push!(B,phyts_b)
write_output(t,phyts_b,dvid_ct,graz_ct,output)
write_output(t,phyts_b,dvid_ct,graz_ct,death_ct,output)
agent_num = size(phyts_b,1)
F = compute_nut_biochem(nutrients, remin)
gtr = compute_source_term(nutrients, velᵇ, g, F)
nutₜ = nut_update(nutrients, consume, g, gtr, ΔT)
write_nut_nc(g, nutₜ, t)
write_nut_cons(g, gtr, nutₜ, velᵇ, t)
write_nut_cons(g, gtr, nutₜ, velᵇ, agent_num, t)
global nutrients = nutₜ;
end

Expand Down Expand Up @@ -148,4 +149,86 @@ open("results/IR.bin", "w") do io
serialize(io, IR)
end

#=t = 1
phyts_a = copy(B[t]) # read data from last time step
phyts_b,dvid_ct,graz_ct,death_ct,consume=phyt_update(t, ΔT, g, phyts_a, nutrients, IR, temp)
velᵇ = read_offline_vels(vfroot,itList,tN,trunc(Int,t*ΔT/3600));
velᵈ = double_grid(velᵇ,g)
agent_move(phyts_a,velᵈ,g,ΔT)
push!(B,phyts_b)
write_output(t,phyts_b,dvid_ct,graz_ct,death_ct,output)
agent_num = size(phyts_b,1)
F = compute_nut_biochem(nutrients, remin)
gtr = compute_source_term(nutrients, velᵇ, g, F)
nutₜ = nut_update(nutrients, consume, g, gtr, ΔT)
global nutrients = nutₜ;
maximum(consume.DON)
findall(isequal(minimum(.POC)),F.POC)
Btest=DataFrame(x=Float64[], y=Float64[], z=Float64[], gen=Int64[], size=Float64[], Cq1=Float64[], Cq2=Float64[], Nq=Float64[], chl=Float64[], sp=Int64[])
Btest1=DataFrame(x=Float64[], y=Float64[], z=Float64[], gen=Int64[], size=Float64[], Cq1=Float64[], Cq2=Float64[], Nq=Float64[], chl=Float64[], sp=Int64[])
Btest2=DataFrame(x=Float64[], y=Float64[], z=Float64[], gen=Int64[], size=Float64[], Cq1=Float64[], Cq2=Float64[], Nq=Float64[], chl=Float64[], sp=Int64[]);
for i in findall(x -> 187 > x >184,phyts_a.x)
push!(Btest,phyts_a[i,:])
end
for i in findall(x -> 262 > x > 260,Btest.y)
push!(Btest1,Btest[i,:])
end
agent_move(Btest1,velᵈ,g,ΔT)
phyts_a = copy(B[7])
Num_phyt = size(phyts_a,1)
chl_num = count_chl(phyts_a, g)
cumsum_chl = cumsum(chl_num, dims = 3);
phyt =phyts_a[615,:]
z = trunc(Int, phyt.z); x = trunc(Int, phyt.x); y = trunc(Int, phyt.y);
DIN = max(0.0, nutrients.DIN[x, y, z])
#compute probabilities of grazing and division
P_graz = rand(Bernoulli(exp(Num_phyt/N*Nsp)*phyt.size/Grz_P))
# Hypothesis: the population of grazers is large enough to graze on phytoplanktons
P_dvi=max(0.0,phyt.size-dvid_size)*1.0e5*rand(Bernoulli(phyt.size/Dvid_P))
PAR = PAR_cal(IR[trunc(Int,t*ΔT/3600)], g.zF[z], cumsum_chl[x, y, z])
PP = PC(PAR,temp[trunc(Int,t*ΔT/3600)],phyt)*phyt.Cq2*ΔT
VN = min(DIN*g.V[x,y,z]/10.0,Nuptake(DIN,phyt)*phyt.Cq2*ΔT)
Dmd_NC = (1+respir_extra(phyt.Cq1))*VN/R_NC
Res2 = k_respir(phyt.size)*phyt.Cq2*ΔT
ρ_chl = chl_sync(phyt,PC(PAR,temp[trunc(Int,t*ΔT/3600)],phyt),IR[trunc(Int,t*ΔT/3600)])
if PP > Dmd_NC
ExuC = (PP - Dmd_NC)*FracExuC
CostC= Dmd_NC
SynC = VN/R_NC
else
CostC= PP
SynC = PP/(1+respir_extra(phyt.Cq1))
ExuC = 0.0
VN = SynC*R_NC
end #exudation
VN
DIN
dCq1 = PP - CostC - ExuC
dCq2 = SynC - Res2
dNq = VN - Res2*R_NC
dsize= dCq2/phyt.Cq2
phyt.Cq2+phyt.Cq1 ≥ Cmin*Nn
consume.DIC[x, y, z] = consume.DIC[x, y, z] + Res2 + CostC - SynC
consume.DIN[x, y, z] = consume.DIN[x, y, z] - VN + Res2*R_NC
consume.DOC[x, y, z] = consume.DOC[x, y, z] + ExuC
consume.DIC[x, y, z]=#


12 changes: 6 additions & 6 deletions nutrient_processes.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,12 +50,12 @@ function nut_update(nutrients, consume, g, gtr, ΔT)
for k in 1:g.Nz
for j in 1:g.Ny
for i in 1:g.Nx
nutₜ.DIC[i, j, k] = nutrients.DIC[i, j, k] + gtr.DIC[i, j, k] * ΔT + consume.DIC[i, j, k]
nutₜ.DIN[i, j, k] = nutrients.DIN[i, j, k] + gtr.DIN[i, j, k] * ΔT + consume.DIN[i, j, k]
nutₜ.DOC[i, j, k] = nutrients.DOC[i, j, k] + gtr.DOC[i, j, k] * ΔT + consume.DOC[i, j, k]
nutₜ.DON[i, j, k] = nutrients.DON[i, j, k] + gtr.DON[i, j, k] * ΔT + consume.DON[i, j, k]
nutₜ.POC[i, j, k] = nutrients.POC[i, j, k] + gtr.POC[i, j, k] * ΔT + consume.POC[i, j, k]
nutₜ.PON[i, j, k] = nutrients.PON[i, j, k] + gtr.PON[i, j, k] * ΔT + consume.PON[i, j, k]
nutₜ.DIC[i, j, k] = nutrients.DIC[i, j, k] + gtr.DIC[i, j, k] * ΔT + consume.DIC[i, j, k]/g.V[i, j, k]
nutₜ.DIN[i, j, k] = nutrients.DIN[i, j, k] + gtr.DIN[i, j, k] * ΔT + consume.DIN[i, j, k]/g.V[i, j, k]
nutₜ.DOC[i, j, k] = nutrients.DOC[i, j, k] + gtr.DOC[i, j, k] * ΔT + consume.DOC[i, j, k]/g.V[i, j, k]
nutₜ.DON[i, j, k] = nutrients.DON[i, j, k] + gtr.DON[i, j, k] * ΔT + consume.DON[i, j, k]/g.V[i, j, k]
nutₜ.POC[i, j, k] = nutrients.POC[i, j, k] + gtr.POC[i, j, k] * ΔT + consume.POC[i, j, k]/g.V[i, j, k]
nutₜ.PON[i, j, k] = nutrients.PON[i, j, k] + gtr.PON[i, j, k] * ΔT + consume.PON[i, j, k]/g.V[i, j, k]
end
end
end
Expand Down
33 changes: 18 additions & 15 deletions parameters.jl
View file
Original file line number Diff line number Diff line change
@@ -1,27 +1,30 @@
####################################### Parameters ########################################
PCmax = [2.7, 2.5] ./86400 # Maximum primary production rate (per second)
Chl2N = 3.0 # Maximum Chla:N ratio in phytoplankton
R_NC = 16/106 # N:C ratio in cell biomass, should be lower than 16:106
Cmin = 7.0e-7 # Minimum C quota (mmolC)
PC_b = [-0.15, -0.2] # Shape parameter
Chl2N = 3.0 # Maximum Chla:N ratio in phytoplankton
R_NC = 16/106 # N:C ratio in cell biomass, should be lower than 16:106
Cmin = 5.0e-13 # Minimum C quota (mmolC)
dvdcount = 0
α = 0.030 # Irradiance absorption coeff (m^2/gChl)
katten_w = 0.046 # PAR attenuation (/m)
katten_c = 0.04 # PAR attenuation (/cell)
Cquota = [1.8e-11, 9.0e-11] # Average C quota in cell (mmolC).
Grz_P = 1500 # phyt.size/Grz_P is the probability to be grazed
dvid_size = 1.6 # relative cell size a cell can start divide
Dvid_P = 6.0 # should be greater than 2.0,phyt.size/Dvid_P is the probability to divide
α = 0.030 # Irradiance absorption coeff (m^2/gChl)
katten_w = 0.046 # PAR attenuation (/m)
katten_c = 0.14 # PAR attenuation (/chl/m)
Cquota = [1.8e-11, 1.8e-10] # Average C quota in cell (mmolC).
Grz_P = 2000 # phyt.size/Grz_P is the probability to be grazed
dvid_size = 1.6 # relative cell size a cell can start divide
Dvid_P = 6.0 # should be greater than 2.0,phyt.size/Dvid_P is the probability to divide

Tempref = 293.15 # reference temperature in K
TempAe = -4000.0 # Arrenhius equation
TempCoeff = 0.8 # Arrenhius equation

VNmax = [1.5, 1.2] ./ 86400 # Maximum N uptake rate (per second)
Nqmax_a = 0.5 # Maximum N quota in cell (mmol N/mmol C)
Nqmin_a = 0.13 # Minimum N quota in cell (mmol N/mmol C)
KsatN = 0.05 # Half-saturation coeff
VNmax = [1.2, 1.0] ./ 86400 # Maximum N uptake rate (per second)
Nqmax = 0.5 # Maximum N quota in cell (mmol N/mmol C)
Nqmin = 0.13 # Minimum N quota in cell (mmol N/mmol C)
KsatN = 0.05 # Half-saturation coeff
VN_b = [-0.15, -0.2] # Shape parameter

respir_a = 1.0e-22 # Respiration ratio
respir_a = 1.0e-20 # Respiration ratio
respir_ex= 1.0e-13 # Extra cost of C for biosynthesis
respir_b = 0.93 # Shape parameter

grazFracC = 0.7 # Fraction goes into dissolved organic pool
Expand Down
Loading

0 comments on commit dd045d3

Please sign in to comment.