Skip to content

sync dev with main #82

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
andybeet wants to merge 7 commits into
base: main
Choose a base branch
from
Open

sync dev with main #82

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
b7b6bc4
Fixed timing issue with contaminant transfer
Sep 6, 2026
671e5f0
Merge branch 'contaminanttiming' of https://github.com/Atlantis-Ecosy…
hmorzaria Sep 6, 2026
2b537e2
eliminated debugging function Debug_Log_Initial_Concentrations
hmorzaria Sep 6, 2026
dc73474
Commented remaining debug statements
hmorzaria Sep 6, 2026
7471bad
Cleaned comments
hmorzaria Sep 6, 2026
a7ebd5d
Merge pull request #79 from Atlantis-Ecosystem-Model/main
eafulton Sep 7, 2026
4e078de
Merge pull request #78 from Atlantis-Ecosystem-Model/contaminant-uptake
raishalovindeer Sep 10, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Binary file added .RData
View file
Open in desktop
Binary file not shown.
136 changes: 136 additions & 0 deletions .Rhistory
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,136 @@
rm(list = ls()) # clear memory
library("devtools")
library("ReactiveAtlantis")
library("proj4")
library("here")
####### Compare outputs and Biomass visualization #######
nc.current <- here::here("example","outputFolder_original","outputSETAS.nc")
nc.old <- here::here("example","outputFolder_mum","outputSETAS.nc")
grp.csv <- 'SETasGroupsDem.csv'
bgm.file <- 'VMPA_setas.bgm'
cum.depths <- c(0,20,50,100,250,700,2000) ## This should be the cummulative depth of your model
## individual file
compare(nc.current, nc.out.old = NULL, grp.csv, bgm.file, cum.depths)
grp.csv <- here::here("example","SETasGroupsDem.csv")
bgm.file <- here::here("example","VMPA_setas.bgm")
cum.depths <- c(0,20,50,100,250,700,2000) ## This should be the cummulative depth of your model
## individual file
compare(nc.current, nc.out.old = NULL, grp.csv, bgm.file, cum.depths)
## compare to previous run
compare(nc.current, nc.old, grp.csv, bgm.file, cum.depths)
biom <- here::here("example","outputFolder_original","outputSETASBiomIndx.txt")
diet.file <- here::here("example","outputFolder_original","outputSETASDietCheck.txt")
bio.age <- here::here("example","outputFolder_original","outputSETASAgeBiomIndx.txt") ## optional file. just if you want to check the predation by age
grp.csv <- here::here("example","SETasGroupsDem.csv")
## Predation by Age
predation(biom, grp.csv, diet.file, bio.age)
## No predation by Age
predation(biom, grp.csv, diet.file, bio.age = NULL)
prm.file <- here::here("example","outputFolder_original","VMPA_setas_biol_fishing_Trunk.prm")
nc.initial <- here::here("example","outputFolder_original","INIT_VMPA_Jan2015.nc")
grp.csv <- here::here("example","SETasGroupsDem.csv")
bgm.file <- here::here("example","VMPA_setas.bgm")
cum.depths <- c(0,20,50,100,250,700,2000) ## This should be the cummulative depth of your model
feeding.mat(prm.file, grp.csv, nc.initial, bgm.file, cum.depths)
prm.file <- here::here("example","VMPA_setas_biol_fishing_Trunk.prm")
nc.initial <- here::here("example","INIT_VMPA_Jan2015.nc")
grp.csv <- here::here("example","SETasGroupsDem.csv")
bgm.file <- here::here("example","VMPA_setas.bgm")
cum.depths <- c(0,20,50,100,250,700,2000) ## This should be the cummulative depth of your model
feeding.mat(prm.file, grp.csv, nc.initial, bgm.file, cum.depths)
grp.csv <- here::here("example","SETasGroupsDem.csv")
prm.file <- here::here("example","VMPA_setas_biol_fishing_Trunk.prm")
diet.file <- here::here("example","outputFolder_original","outputSETASDietCheck.txt")
food.web(diet.file, grp.csv)
## optional you can explore the food web by polygon
food.web(diet.file, grp.file, diet.file.bypol)
food.web(diet.file, grp.csv)
nc.initial <- here::here("example","INIT_VMPA_Jan2015.nc")
nc.current <- here::here("example","outputFolder_original","outputSETAS.nc")
grp.csv <- here::here("example","SETasGroupsDem.csv")
prm.file <- here::here("example","VMPA_setas_biol_fishing_Trunk.prm")
growth.pp(nc.initial, grp.csv, prm.file, nc.current)
rm(list = ls()) # clear memory
library("devtools")
library("ReactiveAtlantis")
library("proj4")
library("here")
####### Compare outputs and Biomass visualization #######
nc.current <- here::here("example","outputFolder_original","outputSETAS.nc")
nc.old <- here::here("example","outputFolder_mum","outputSETAS.nc")
grp.csv <- here::here("example","SETasGroupsDem.csv")
bgm.file <- here::here("example","VMPA_setas.bgm")
cum.depths <- c(0,20,50,100,250,700,2000) ## This should be the cummulative depth of your model
## individual file
compare(nc.current, nc.out.old = NULL, grp.csv, bgm.file, cum.depths)
## compare to previous run
compare(nc.current, nc.old, grp.csv, bgm.file, cum.depths)
####### Predation analysis from the Atlantis output #######
biom <- here::here("example","outputFolder_original","outputSETASBiomIndx.txt")
diet.file <- here::here("example","outputFolder_original","outputSETASDietCheck.txt")
bio.age <- here::here("example","outputFolder_original","outputSETASAgeBiomIndx.txt") ## optional file. just if you want to check the predation by age
grp.csv <- here::here("example","SETasGroupsDem.csv")
## Predation by Age
predation(biom, grp.csv, diet.file, bio.age)
## No predation by Age
predation(biom, grp.csv, diet.file, bio.age = NULL)
prm.file <- here::here("example","VMPA_setas_biol_fishing_Trunk.prm")
nc.initial <- here::here("example","INIT_VMPA_Jan2015.nc")
grp.csv <- here::here("example","SETasGroupsDem.csv")
bgm.file <- here::here("example","VMPA_setas.bgm")
cum.depths <- c(0,20,50,100,250,700,2000) ## This should be the cummulative depth of your model
feeding.mat(prm.file, grp.csv, nc.initial, bgm.file, cum.depths)
feeding.mat(prm.file, grp.csv, nc.initial, bgm.file, cum.depths)
grp.csv <- here::here("example","SETasGroupsDem.csv")
prm.file <- here::here("example","VMPA_setas_biol_fishing_Trunk.prm")
diet.file <- here::here("example","outputFolder_original","outputSETASDietCheck.txt")
food.web(diet.file, grp.csv)
nc.initial <- here::here("example","INIT_VMPA_Jan2015.nc")
nc.current <- here::here("example","outputFolder_original","outputSETAS.nc")
grp.csv <- here::here("example","SETasGroupsDem.csv")
prm.file <- here::here("example","VMPA_setas_biol_fishing_Trunk.prm")
growth.pp(nc.initial, grp.csv, prm.file, nc.current)
#
nc.initial <- here::here("example","INIT_VMPA_Jan2015.nc")
nc.current <- here::here("example","outputFolder_original","outputSETAS.nc")
yoy.file <- here::here("example","outputFolder_original","outputSETASYOY.txt")
grp.csv <- here::here("example","SETasGroupsDem.csv")
prm.file <- here::here("example","VMPA_setas_biol_fishing_Trunk.prm")
recruitment.cal(nc.initial, nc.current, yoy.file, grp.csv, prm.file)
rm(list = ls()) # clear memory
rm(list = ls()) # clear memory
library("devtools")
library("ReactiveAtlantis")
library("proj4")
library("here")
####### Compare outputs and Biomass visualization #######
nc.current <- here::here("example","outputFolder_original","outputSETAS.nc")
nc.old <- here::here("example","outputFolder_mum","outputSETAS.nc")
grp.csv <- here::here("example","SETasGroupsDem.csv")
bgm.file <- here::here("example","VMPA_setas.bgm")
cum.depths <- c(0,20,50,100,250,700,2000) ## This should be the cummulative depth of your model
## individual file
compare(nc.current, nc.out.old = NULL, grp.csv, bgm.file, cum.depths)
compare(nc.current, nc.old, grp.csv, bgm.file, cum.depths)
biom <- here::here("example","outputFolder_original","outputSETASBiomIndx.txt")
diet.file <- here::here("example","outputFolder_original","outputSETASDietCheck.txt")
bio.age <- here::here("example","outputFolder_original","outputSETASAgeBiomIndx.txt") ## optional file. just if you want to check the predation by age
grp.csv <- here::here("example","SETasGroupsDem.csv")
## Predation by Age
predation(biom, grp.csv, diet.file, bio.age)
####### Atlantis food web and trophic level composition #######
grp.csv <- here::here("example","SETasGroupsDem.csv")
prm.file <- here::here("example","VMPA_setas_biol_fishing_Trunk.prm")
diet.file <- here::here("example","outputFolder_original","outputSETASDietCheck.txt")
food.web(diet.file, grp.csv)
biom <- here::here("example","outputFolder_original","outputSETASBiomIndx.txt")
diet.file <- here::here("example","outputFolder_original","outputSETASDietCheck.txt")
bio.age <- here::here("example","outputFolder_original","outputSETASAgeBiomIndx.txt") ## optional file. just if you want to check the predation by age
grp.csv <- here::here("example","SETasGroupsDem.csv")
## Predation by Age
predation(biom, grp.csv, diet.file, bio.age)
claude
1850/365
source("~/psatlantismodel/PS_runs_PCBs.R")
source("~/psatlantismodel/PS_runs_PCBs.R")
source("~/psatlantismodel/PS_runs_PCBs.R")
99 changes: 70 additions & 29 deletions atlantis/atecology/atContaminants.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,12 +31,8 @@
#include "atecology.h"

static FILE *contaminantContactFile;
/**
* Free up the contaminant structure.
*
*
*
*/


void Free_Contaminants(MSEBoxModel *bm) {
int cIndex;

Expand Down Expand Up @@ -422,10 +418,15 @@ void Move_Vert_Contaminated(MSEBoxModel *bm, int sp, int cohort, double ****this
*/

void Change_Contaminant_Levels(MSEBoxModel *bm, BoxLayerValues *boxLayerInfo, HABITAT_TYPES habitat, double dtsz, int offset) {
int cIndex;
int cIndex, sp, cohort;
double cLevel;
double *tracerArray = getTracerArray(boxLayerInfo, habitat);

/* DEBUG: Log when contaminant processing happens
if (bm->dayt < 2.0 && habitat == WC) {
fprintf(bm->logFile, "DEBUG_TIMING: Change_Contaminant_Levels START for habitat WC at time=%e\n", bm->dayt);
} */

for (cIndex = 0; cIndex < bm->num_contaminants; cIndex++) {
/* Grab the level in the water column or the sediment */
cLevel = tracerArray[bm->contaminantStructure[cIndex]->contaminant_tracer];
Expand All @@ -436,6 +437,11 @@ void Change_Contaminant_Levels(MSEBoxModel *bm, BoxLayerValues *boxLayerInfo, HA
}
Calculate_Contaminants_Flux(bm, boxLayerInfo, habitat);

/* DEBUG: Log when contaminant processing ends
if (bm->dayt < 2.0 && habitat == WC) {
fprintf(bm->logFile, "DEBUG_TIMING: Change_Contaminant_Levels END for habitat WC at time=%e\n", bm->dayt);
} */

return;
}

Expand Down Expand Up @@ -742,32 +748,46 @@ int Species_Contaminant_Uptake(MSEBoxModel *bm, BoxLayerValues *boxLayerInfo, HA

int sp, pid, cohort = 0;
double cLevel = 0, uptake_rate = 0, cUptake = 0, cGroupLevel = 0, Cnew = 0, propLevel = 0, avoidLevel = 0, chronicLevel = 0, testLevel = 0, propContam = 0, time_step = 0;
double *tracerArray;
double *environmentTracerArray, *organismTracerArray;
CONTAMINANT_UPTAKE_OPTION uptake_option;

/* For EPIFAUNA groups: environmental concentration comes from WC, but organism concentration lives in EPIFAUNA */
if (habitat == EPIFAUNA) {
tracerArray = getTracerArray(boxLayerInfo, WC);
environmentTracerArray = getTracerArray(boxLayerInfo, WC); /* Water column concentration (forced from outside) */
organismTracerArray = getTracerArray(boxLayerInfo, EPIFAUNA); /* Organism concentration in EPIFAUNA */
} else {
tracerArray = getTracerArray(boxLayerInfo, habitat);
/* For other habitats, both come from the same habitat */
environmentTracerArray = getTracerArray(boxLayerInfo, habitat);
organismTracerArray = getTracerArray(boxLayerInfo, habitat);
}


/* For each contaminant calculate the uptake */

/* Grab the level in the water column or the sediment */
cLevel = tracerArray[bm->contaminantStructure[cIndex]->contaminant_tracer];
/* Grab the level in the water column or the sediment (environmental concentration) */
cLevel = environmentTracerArray[bm->contaminantStructure[cIndex]->contaminant_tracer];

if(cLevel > bm->min_pool_cont){
//fprintf(bm->logFile, "time %e, box %d, layer %d, cLevel = %e\n", bm->dayt, bm->current_box, bm->current_layer, cLevel);
for (sp = 0; sp < bm->K_num_tot_sp; sp++) {

/* Contact can occur for all active groups that are present in this habitat */
if ((FunctGroupArray[sp].speciesParams[flag_id] == TRUE) && (FunctGroupArray[sp].isDetritus == FALSE) && (FunctGroupArray[sp].habitatCoeffs[habitat] > 0)) {

uptake_option = (CONTAMINANT_UPTAKE_OPTION)bm->contaminantStructure[cIndex]->sp_uptake_option[sp];
uptake_rate = bm->contaminantStructure[cIndex]->sp_uptake_rate[sp];

for(cohort = 0; cohort < FunctGroupArray[sp].numCohortsXnumGenes; cohort++){
/* The current concentration in the group */
cGroupLevel = tracerArray[FunctGroupArray[sp].contaminantTracers[cohort][cIndex]];
int h;
/* The current concentration in the group (from organism's habitat) */
cGroupLevel = organismTracerArray[FunctGroupArray[sp].contaminantTracers[cohort][cIndex]];
/* Add any transfers from this timestep that haven't been integrated into tracer yet */
/* sp_transfer is indexed by prey habitat, not predator habitat, so sum across all habitats */
for (h = 0; h < bm->num_active_habitats; h++) {
if (bm->contaminantStructure[cIndex]->sp_transfer[sp][cohort][h] > 0) {
cGroupLevel += bm->contaminantStructure[cIndex]->sp_transfer[sp][cohort][h];
}
}
cUptake = 0;

switch (uptake_option) {
Expand All @@ -791,12 +811,12 @@ int Species_Contaminant_Uptake(MSEBoxModel *bm, BoxLayerValues *boxLayerInfo, HA
}

if ((time_step > 0) && (uptake_rate > 0)) {
/* Have the calculate the actual new tracer value and from this we calculate the flux values
/* Have the calculate the actual new tracer value and from this we calculate the flux values
contam_sig_uptake_const set to 0.99 by default
*/

//fprintf(bm->logFile,"Time: %e box %d-%d contam_sig_uptake_const: %e\n", bm->dayt, bm->current_box, bm->current_layer, bm->contam_sig_uptake_const);

Cnew = pow((pow(cLevel,(1.0 - bm->contam_sig_uptake_const))- (pow(cLevel,(1.0 - bm->contam_sig_uptake_const)) - pow(cGroupLevel, (1.0 - bm->contam_sig_uptake_const))) * exp(-uptake_rate * time_step * (bm->contam_sig_uptake_const + 1.0))),(1.0 / (1.0 - bm->contam_sig_uptake_const)));
cUptake = (Cnew - cGroupLevel) / time_step;
//cUptake = 0;
Expand Down Expand Up @@ -830,9 +850,9 @@ int Species_Contaminant_Uptake(MSEBoxModel *bm, BoxLayerValues *boxLayerInfo, HA
break;
}


//fprintf(bm->logFile,"Time: %e box%d-%d %s-%d vs %s has uptake_option: %d cUptake: %e uptake_rate: %e cLevel: %e\n", bm->dayt, bm->current_box, bm->current_layer, FunctGroupArray[sp].groupCode, cohort, bm->contaminantStructure[cIndex]->contaminant_name, uptake_option, cUptake, uptake_rate, cLevel);

bm->contaminantStructure[cIndex]->sp_uptake[sp][cohort][habitat] = cUptake;

// Initiate exposure time
Expand Down Expand Up @@ -940,7 +960,7 @@ int Init_Contaminant_Transfer_Values(MSEBoxModel *bm) {
int Group_Transfer_Contaminant(MSEBoxModel *bm, BoxLayerValues *boxLayerInfo, HABITAT_TYPES globalHabitat, HABITAT_TYPES habitat, int toGuild, int toCohort, int fromGuild, int fromCohort, double amountTransfer, double ***spSPinfo, double initialBiomass, double dtsz, int need_prop, int caseGTC) {

int cIndex, pid, this_habitat;
double *tracerArray;
double *preyTracerArray, *predatorTracerArray;
double cGroupLevel = 0, transfer = 0, totalBiomass, amt_exchanged, toGuild_totalBiomass = 0.0, prop_exchanged, propContam, min_num, this_num;
int isGlobal = (FunctGroupArray[toGuild].diagTol == 2 && it_count == 1);

Expand All @@ -949,15 +969,22 @@ int Group_Transfer_Contaminant(MSEBoxModel *bm, BoxLayerValues *boxLayerInfo, HA
*/
if (amountTransfer == 0.0)
return 0;

if(isnan(dtsz)){
quit("Group_Transfer_Contaminant - from group %s-%d, to group %s-%d, dtsz: %e\n", FunctGroupArray[fromGuild].groupCode, fromCohort, FunctGroupArray[toGuild].groupCode, toCohort, dtsz);
}

/* Get tracer arrays - separate for prey and predator since they may be in different habitats */
if(FunctGroupArray[fromGuild].habitatType == EPIFAUNA){
tracerArray = getTracerArray(boxLayerInfo, EPIFAUNA);
preyTracerArray = getTracerArray(boxLayerInfo, EPIFAUNA);
} else {
preyTracerArray = getTracerArray(boxLayerInfo, habitat);
}

if(FunctGroupArray[toGuild].habitatType == EPIFAUNA){
predatorTracerArray = getTracerArray(boxLayerInfo, EPIFAUNA);
} else {
tracerArray = getTracerArray(boxLayerInfo, habitat);
predatorTracerArray = getTracerArray(boxLayerInfo, habitat);
}

/* Work out the proportion transfered */
Expand Down Expand Up @@ -994,8 +1021,8 @@ int Group_Transfer_Contaminant(MSEBoxModel *bm, BoxLayerValues *boxLayerInfo, HA
quit("Time: %e %s-%d returning nan biomass in box%d-%d for Group_Transfer_Contaminant - SN: %e RN: %e DEN: %e cell_vol: %e\n", bm->dayt, FunctGroupArray[toGuild].groupCode, toCohort, bm->current_box, bm->current_layer, spSPinfo[toGuild][toCohort][SN_id], spSPinfo[toGuild][toCohort][RN_id], spSPinfo[toGuild][toCohort][DEN_id], bm->cell_vol);
}
} else {
toGuild_totalBiomass = tracerArray[FunctGroupArray[toGuild].totNTracers[toCohort]];
toGuild_totalBiomass = predatorTracerArray[FunctGroupArray[toGuild].totNTracers[toCohort]];

if(isnan(toGuild_totalBiomass)) {
quit("Time: %e %s-%d returning nan biomass in box%d-%d for Group_Transfer_Contaminant\n", bm->dayt, FunctGroupArray[toGuild].groupCode, toCohort, bm->current_box, bm->current_layer);
}
Expand All @@ -1004,9 +1031,10 @@ int Group_Transfer_Contaminant(MSEBoxModel *bm, BoxLayerValues *boxLayerInfo, HA


for (cIndex = 0; cIndex < bm->num_contaminants; cIndex++) {
/* The current concentration in the group */
cGroupLevel = tracerArray[FunctGroupArray[fromGuild].contaminantTracers[fromCohort][cIndex]];
/* The current concentration in the prey */
cGroupLevel = preyTracerArray[FunctGroupArray[fromGuild].contaminantTracers[fromCohort][cIndex]];

/* DEBUG: Log feeding with timestamps */
/* This used to have an if-else statement using the test (cGroupLevel > bm->min_pool) but removed so got sediment transferal and accumualtion even if low amounts
Use of Contamflag means it won't get reset in Integrate_Tracer_Variables() */

Expand All @@ -1032,6 +1060,19 @@ int Group_Transfer_Contaminant(MSEBoxModel *bm, BoxLayerValues *boxLayerInfo, HA
bm->contaminantStructure[cIndex]->sp_transfer[toGuild][toCohort][habitat] += transfer;
bm->contaminantStructure[cIndex]->sp_transfer[fromGuild][fromCohort][habitat] -= transfer;

/* CRITICAL FIX: Apply transfer immediately to tracer arrays so predators see updated concentrations */
/* This ensures inverts see prey contaminants immediately after eating, not one timestep later */
/* Validate tracer indices before accessing to prevent out-of-bounds access */
int toGuild_tracer_idx = FunctGroupArray[toGuild].contaminantTracers[toCohort][cIndex];
int fromGuild_tracer_idx = FunctGroupArray[fromGuild].contaminantTracers[fromCohort][cIndex];

if (toGuild_tracer_idx >= 0 && toGuild_tracer_idx < bm->ntracer) {
predatorTracerArray[toGuild_tracer_idx] += transfer;
}
if (fromGuild_tracer_idx >= 0 && fromGuild_tracer_idx < bm->ntracer) {
preyTracerArray[fromGuild_tracer_idx] -= transfer;
}

/**
//if((((toGuild == 54) || (fromGuild == 54)) && (bm->contaminantStructure[cIndex]->sp_transfer[toGuild][toCohort][habitat] > 0.0)) && (cIndex == 3)) {
fprintf(bm->logFile, "prey = %s, to %s-%d gaining %e, cGroupLevel= %e, propTransfer= %e, totalTransfer = %e amountTransfer = %e, totalBiomass= %e\n", FunctGroupArray[fromGuild].groupCode, FunctGroupArray[toGuild].groupCode, toCohort, transfer, cGroupLevel, propTransfer, bm->contaminantStructure[cIndex]->sp_transfer[toGuild][toCohort][habitat], amountTransfer, totalBiomass);
Expand All @@ -1055,7 +1096,7 @@ int Group_Transfer_Contaminant(MSEBoxModel *bm, BoxLayerValues *boxLayerInfo, HA

min_num = 0.0;
if(FunctGroupArray[toGuild].groupAgeType == AGE_STRUCTURED ){
this_num = tracerArray[FunctGroupArray[toGuild].NumsTracers[toCohort]];
this_num = predatorTracerArray[FunctGroupArray[toGuild].NumsTracers[toCohort]];
min_num = 1.0 / (this_num + small_num);
} else {
min_num = bm->min_pool_cont; // Just using this as a proxy for a small number here
Expand Down
Loading
Loading