From dd045d3a18b48760e2415f4dd4455a59bfe1e120 Mon Sep 17 00:00:00 2001 From: zhenwu0728 Date: Sat, 1 Jun 2019 17:26:19 -0400 Subject: [PATCH] physiologycal processes improved --- agent_div.jl | 2 +- model_setup.jl | 24 +++---- model_update.jl | 99 +++++++++++++++++++++++++--- nutrient_processes.jl | 12 ++-- parameters.jl | 33 +++++----- phyt_process.jl | 147 ++++++++++++++++++++++-------------------- utils.jl | 12 ++-- 7 files changed, 212 insertions(+), 117 deletions(-) diff --git a/agent_div.jl b/agent_div.jl index ab22a87f..1622d19d 100644 --- a/agent_div.jl +++ b/agent_div.jl @@ -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 diff --git a/model_setup.jl b/model_setup.jl index def2e4b0..be1bf6a7 100644 --- a/model_setup.jl +++ b/model_setup.jl @@ -3,8 +3,8 @@ 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))) @@ -12,15 +12,15 @@ function setup_agents(N::Int64,Cquota::Array,Nn::Int64,mean::Float64,var::Float6 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))) @@ -28,8 +28,8 @@ function setup_agents(N::Int64,Cquota::Array,Nn::Int64,mean::Float64,var::Float6 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 @@ -46,9 +46,9 @@ 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] @@ -56,7 +56,7 @@ function setup_nutrients(g,nut) 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] diff --git a/model_update.jl b/model_update.jl index 23afb3a1..f6af1240 100644 --- a/model_update.jl +++ b/model_update.jl @@ -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 @@ -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]=# + diff --git a/nutrient_processes.jl b/nutrient_processes.jl index a36f2966..d23fb9a5 100644 --- a/nutrient_processes.jl +++ b/nutrient_processes.jl @@ -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 diff --git a/parameters.jl b/parameters.jl index da5ef333..627cf122 100644 --- a/parameters.jl +++ b/parameters.jl @@ -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 diff --git a/phyt_process.jl b/phyt_process.jl index 287a7622..d0e46d56 100644 --- a/phyt_process.jl +++ b/phyt_process.jl @@ -1,8 +1,8 @@ ########################################### # define functions of pysiology processes # ########################################### -k_respir(Cq2) = respir_a*(12.0e9*Cq2)^respir_b/Cq2 -respir(Cq2) = respir_a*(12.0e9*Cq2)^respir_b +respir_extra(Cq1) = respir_ex*(12.0e9*Cq1/Nn)^respir_b # extra cost for biosynthesis, related to Cq1 +k_respir(size) = respir_a*(size)^respir_b # respiration rate, unit: per second function daynight(t::Int64, IR) if IR[trunc(Int,t*ΔT/3600)] < 5.0 @@ -12,8 +12,8 @@ function daynight(t::Int64, IR) end end -function PAR_cal(I, z, cumsum_cell) - atten = (katten_w + katten_c * cumsum_cell)*(-z) +function PAR_cal(I, z, cumsum_chl) + atten = (katten_w + katten_c * cumsum_chl)*(-z) PAR = α*I*(1.0 - exp(-atten))/atten return PAR end @@ -21,16 +21,17 @@ end function PC(PAR, Temp, phyt) Tempstd = exp(TempAe*(1.0/(Temp+273.15)-1.0/Tempref)) photoTempFunc = TempCoeff*max(1.0e-10,Tempstd) - PC = PCmax[phyt.sp]*photoTempFunc*(1-exp(-PAR*phyt.chl/(phyt.Cq2*PCmax[phyt.sp]*photoTempFunc))) + PCm = PCmax[phyt.sp]*phyt.size^PC_b[phyt.sp] + PC = PCm*photoTempFunc*(1-exp(-PAR*phyt.chl/(phyt.Cq2*PCmax[phyt.sp]*photoTempFunc))) return PC end function Nuptake(Nit, phyt) - Nqmax = Nqmax_a*phyt.size*Cquota[phyt.sp]*Nn - Nqmin = Nqmin_a*phyt.size*Cquota[phyt.sp]*Nn + Qn = phyt.Nq/(phyt.Cq1+phyt.Cq2) #In-Cell N uptake limitation - regQ = max(0.0,min(1.0,(Nqmax-phyt.Nq)/(Nqmax-Nqmin))) - Nuptake = VNmax[phyt.sp]*Nit/(Nit+KsatN)*regQ + regQ = max(0.0,min(1.0,(Nqmax-Qn)/(Nqmax-Nqmin))) + VNm = VNmax[phyt.sp]*phyt.size^VN_b[phyt.sp] + Nuptake = VNm*Nit/(Nit+KsatN)*regQ return Nuptake end @@ -38,7 +39,7 @@ function chl_sync(phyt,PP,I) if I > 0 ρ_chl = PP/(α*I*phyt.chl/phyt.Cq2) else - ρ_chl = 1.0 + ρ_chl = 0.0 end return ρ_chl end @@ -77,8 +78,8 @@ function phyt_update(t::Int64, ΔT::Int64, g, phyts_a, nutrients, IR, temp) # load nutrients dvid_ct = 0; graz_ct = 0; death_ct = 0 Num_phyt = size(phyts_a,1) - cell_num = count_chl(phyts_a, g) - cumsum_cell = cumsum(cell_num, dims = 3) + chl_num = count_chl(phyts_a, g) + cumsum_chl = cumsum(chl_num, dims = 3) #set up a dataframe to record all updated agents phyts_b = DataFrame(x=Float64[], y=Float64[], z=Float64[], gen=Int64[], size=Float64[], Cq1=Float64[], Cq2=Float64[], Nq=Float64[], chl=Float64[], sp=Int64[]) # @@ -93,72 +94,80 @@ function phyt_update(t::Int64, ΔT::Int64, g, phyts_a, nutrients, IR, temp) 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_cell[x, y, z]) - PP = PC(PAR,temp[trunc(Int,t*ΔT/3600)],phyt)*Cquota[phyt.sp]*phyt.size*ΔT*Nn - VN = min(DIN*g.V[x,y,z]/10.0,Nuptake(DIN,phyt)*Cquota[phyt.sp]*phyt.size*ΔT*Nn) - Dmd_NC = (1+k_respir(phyt.Cq1))*VN/R_NC - Res2 = respir(phyt.Cq2)*ΔT + 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 P_graz < 1 #not grazed - if phyt.Cq2+phyt.Cq1 ≤ Cmin # natural death - consume.DOC[x, y, z] = consume.DOC[x, y, z] + (phyt.Cq1+phyt.Cq2)*mortFracC - consume.DON[x, y, z] = consume.DON[x, y, z] + phyt.Nq*mortFracN - consume.POC[x, y, z] = consume.POC[x, y, z] + (phyt.Cq1+phyt.Cq2)*(1.0 - mortFracC) - consume.PON[x, y, z] = consume.PON[x, y, z] + phyt.Nq*(1.0 - mortFracN) - death_ct += 1 - else # not natural death - if daynight(t,IR) # day - if PP > Dmd_NC - ExuC = (PP - Dmd_NC)*FracExuC - CostC= Dmd_NC - SynC = VN/R_NC - else - CostC= PP - SynC = PP/(1+k_respir(phyt.Cq1)) - ExuC = 0.0 - VN = SynC*R_NC - end #exudation + if P_graz == false #not grazed + if daynight(t,IR) # day + 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 + dCq1 = PP - CostC - ExuC + dCq2 = SynC - Res2 + dNq = VN - Res2*R_NC + dsize= dCq2/phyt.Cq2 + phyt.Cq1 = max(Cmin*Nn/10.0, phyt.Cq1 + dCq1) + phyt.Cq2 = max(Cmin*Nn/10.0, phyt.Cq2 + dCq2) + phyt.Nq = max(Cmin*Nn*R_NC/10.0, phyt.Nq + dNq) + phyt.size= phyt.size+dsize + phyt.chl = phyt.chl + ρ_chl*VN*Chl2N + if phyt.Cq2+phyt.Cq1 ≥ Cmin*Nn # not natural death + push!(phyts_b,phyt) + else # natural death + consume.DOC[x, y, z] = consume.DOC[x, y, z] + (phyt.Cq1+phyt.Cq2)*mortFracC + consume.DON[x, y, z] = consume.DON[x, y, z] + phyt.Nq*mortFracN + consume.POC[x, y, z] = consume.POC[x, y, z] + (phyt.Cq1+phyt.Cq2)*(1.0 - mortFracC) + consume.PON[x, y, z] = consume.PON[x, y, z] + phyt.Nq*(1.0 - mortFracN) + death_ct += 1 + end # naturan death + 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 + else # night + if P_dvi < 1 # not divide + CostC= min(0.2*phyt.Cq1,Dmd_NC) + SynC = min(VN/R_NC,0.2*phyt.Cq1/(1+respir_extra(phyt.Cq1))) + ExuC = 0.0 + VN = SynC*R_NC dCq1 = PP - CostC - ExuC dCq2 = SynC - Res2 dNq = VN - Res2*R_NC dsize= dCq2/phyt.Cq2 - phyt.Cq1 = max(0.0,phyt.Cq1 + dCq1) - phyt.Cq2 = max(0.0,phyt.Cq2 + dCq2) - phyt.Nq = max(0.0,phyt.Nq + dNq) + phyt.Cq1 = max(Cmin*Nn/10.0,phyt.Cq1 + dCq1) + phyt.Cq2 = max(Cmin*Nn/10.0,phyt.Cq2 + dCq2) + phyt.Nq = max(Cmin*Nn/10.0,phyt.Nq + dNq) phyt.size= phyt.size+dsize phyt.chl = phyt.chl + ρ_chl*VN*Chl2N - push!(phyts_b,phyt) + if phyt.Cq2+phyt.Cq1 ≥ Cmin*Nn # not natural death + push!(phyts_b,phyt) + else # natural death + consume.DOC[x, y, z] = consume.DOC[x, y, z] + (phyt.Cq1+phyt.Cq2)*mortFracC + consume.DON[x, y, z] = consume.DON[x, y, z] + phyt.Nq*mortFracN + consume.POC[x, y, z] = consume.POC[x, y, z] + (phyt.Cq1+phyt.Cq2)*(1.0 - mortFracC) + consume.PON[x, y, z] = consume.PON[x, y, z] + phyt.Nq*(1.0 - mortFracN) + death_ct += 1 + end # natural death 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 - else # night - if P_dvi < 1 # not divide - CostC= min(0.4*phyt.Cq1,Dmd_NC) - SynC = min(VN/R_NC,0.4*phyt.Cq1/(1+k_respir(phyt.Cq1))) - ExuC = 0.0 - VN = SynC*R_NC - dCq1 = PP - CostC - ExuC - dCq2 = SynC - Res2 - dNq = VN - Res2*R_NC - dsize= dCq2/phyt.Cq2 - phyt.Cq1 = max(0.0,phyt.Cq1 + dCq1) - phyt.Cq2 = max(0.0,phyt.Cq2 + dCq2) - phyt.Nq = max(0.0,phyt.Nq + dNq) - phyt.size= phyt.size+dsize - phyt.chl = phyt.chl + ρ_chl*VN*Chl2N - push!(phyts_b,phyt) - 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 - else #divide - dvid_ct += 2 - global dvdcount += 2 - phyts2 = divide(phyt) - append!(phyts_b,phyts2) - consume.DIC[x, y, z] = consume.DIC[x, y, z] + phyt.Cq2*0.1 # consume C when cell is divided - end # divide - end # day night? - end # natural death + else #divide + dvid_ct += 2 + global dvdcount += 2 + phyts2 = divide(phyt) + append!(phyts_b,phyts2) + consume.DIC[x, y, z] = consume.DIC[x, y, z] + phyt.Cq2*0.1 # consume C when cell is divided + end # divide + end # day night? else #grazed graz_ct += 1 consume.DOC[x, y, z] = consume.DOC[x, y, z] + (phyt.Cq1+phyt.Cq2)*grazFracC*0.5 @@ -167,5 +176,5 @@ function phyt_update(t::Int64, ΔT::Int64, g, phyts_a, nutrients, IR, temp) consume.PON[x, y, z] = consume.PON[x, y, z] + phyt.Nq*(1.0 - grazFracN)*0.5 end # graze end # while loop to traverse the array of agents - return phyts_b,dvid_ct,graz_ct,consume + return phyts_b,dvid_ct,graz_ct,death_ct,consume end # for loop of time diff --git a/utils.jl b/utils.jl index 783976f5..169742cd 100644 --- a/utils.jl +++ b/utils.jl @@ -96,11 +96,11 @@ function grid_offline(fieldroot::String) end function create_output(B::Array{DataFrame,1}) - output = DataFrame(time=0, gen_ave=mean(B[1].gen), spec_ave = mean(B[1].sp), Cq1_ave=mean(B[1].Cq1), Cq2_ave=mean(B[1].Cq2), Nq_ave=mean(B[1].Nq), size_ave=mean(B[1].size), chl_ave=mean(B[1].chl), Population=size(B[1],1), dvid=0, graz=0) + output = DataFrame(time=0, gen_ave=mean(B[1].gen), spec_ave = mean(B[1].sp), Cq1_ave=mean(B[1].Cq1), Cq2_ave=mean(B[1].Cq2), Nq_ave=mean(B[1].Nq), size_ave=mean(B[1].size), chl_ave=mean(B[1].chl), Population=size(B[1],1), dvid=0, graz=0,death=0) return output end -function write_output(t,phyts_b,dvid_ct,graz_ct,output) +function write_output(t,phyts_b,dvid_ct,graz_ct,death_ct,output) # summary of current step gen_ave=mean(phyts_b.gen) spec_ave=mean(phyts_b.sp) @@ -109,7 +109,7 @@ function write_output(t,phyts_b,dvid_ct,graz_ct,output) Nq_ave=mean(phyts_b.Nq) size_ave=mean(phyts_b.size) chl_ave=mean(phyts_b.chl) - push!(output,(time=t, gen_ave=gen_ave, spec_ave=spec_ave, Cq1_ave=Cq1_ave, Cq2_ave=Cq2_ave, Nq_ave=Nq_ave, size_ave=size_ave, chl_ave=chl_ave, Population=size(phyts_b,1), dvid=dvid_ct, graz=graz_ct)) + push!(output,(time=t, gen_ave=gen_ave, spec_ave=spec_ave, Cq1_ave=Cq1_ave, Cq2_ave=Cq2_ave, Nq_ave=Nq_ave, size_ave=size_ave, chl_ave=chl_ave, Population=size(phyts_b,1), dvid=dvid_ct, graz=graz_ct, death=death_ct)) return output end @@ -212,7 +212,7 @@ function write_nut_nc(g::grids, nut::nutrient_fields, t::Int64) ncclose(filepath) return nothing end -function write_nut_cons(g::grids, gtr::nutrient_fields, nutₜ::nutrient_fields, vel::velocity, t::Int64) +function write_nut_cons(g::grids, gtr::nutrient_fields, nutₜ::nutrient_fields, vel::velocity, agent_num::Int64, t::Int64) Σgtrⁿ = sum(gtr.DIN .* g.V)+sum(gtr.DON .* g.V)+sum(gtr.PON .* g.V) Σgtrᶜ = sum(gtr.DIC .* g.V)+sum(gtr.DOC .* g.V)+sum(gtr.POC .* g.V) ΣsurFⁿ= sum((nutₜ.DIN[:,:,1]+nutₜ.DON[:,:,1]+nutₜ.PON[:,:,1]) .* g.Az .* vel.w[:,:,1]) @@ -222,6 +222,6 @@ function write_nut_cons(g::grids, gtr::nutrient_fields, nutₜ::nutrient_fields, DINio = open("results/cons_DIN.txt","a"); println(Cio,@sprintf("%3.0f %.16E %.16E %.8E",t,Σgtrᶜ,ΣsurFᶜ,Σgtrᶜ+ΣsurFᶜ)) println(Nio,@sprintf("%3.0f %.16E %.16E %.8E",t,Σgtrⁿ,ΣsurFⁿ,Σgtrⁿ+ΣsurFⁿ)) - println(DINio,@sprintf("%3.0f %.16E",t,ΣDIN)) - close(Cio);close(Nio); + println(DINio,@sprintf("%3.0f %.16E %7.0f",t,ΣDIN,agent_num)) + close(Cio);close(Nio);close(DINio); end