sync dev with main

README.md

@@ -41,7 +41,7 @@ The Atlantis version represents smaller fixes and new features and follows the h
 | Icelantis | Iceland | Iceland | [Jacob Kasper](https://github.com/jacobkasper) | <https://github.com/jacobkasper/Icelantis> |
 |  | South Australia | Australia | [Beth Fulton](https://github.com/eafulton) | NA |
 |  SSAM | Salish Sea | Canada | [Raisha Lovindeer](https://github.com/raishalovindeer) | <https://github.com/SS-Atlantis/salish-sea-atlantis-model> |
-| AMPS | Puget Sound | USA; Canada | [Hem Nalini Morzaria-Luna](https://github.com/hmorzaria) | NA |
+| AMPS | Puget Sound | USA; Canada | [Hem Nalini Morzaria-Luna](https://github.com/hmorzaria) | <https://github.com/atlantis-amps/psatlantismodel> |
 |  | Southern Benguela | South Africa | [Kelly Cisneros](https://github.com/kortegac) | NA |
 |  | Tasman and Golden Bays | New Zealand |  | NA |
 |  | Korean Peninsula | South Korea | [Junggil Seo]() | NA |

atlantis/atecology/atContaminants.c

@@ -31,12 +31,8 @@
 #include "atecology.h"
 
 static FILE *contaminantContactFile;
-/**
- * Free up the contaminant structure.
- *
- *
- *
- */
+
+
 void Free_Contaminants(MSEBoxModel *bm) {
 	int cIndex;
 
@@ -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];
@@ -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;
 }
 
@@ -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) {
@@ -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;
@@ -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
@@ -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);
 
@@ -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 */
@@ -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);
             }
@@ -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() */
 
@@ -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);
@@ -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

atlantis/atecology/atGroupProcesses.c

@@ -158,6 +158,10 @@ static void Update_Detritus(MSEBoxModel *bm, BoxLayerValues *boxLayerInfo, HABIT
             //fprintf(bm->logFile,"Calling Group_Transfer_Contaminant from Update_Detritus WC WC - %s-%d with initialBiomass: %e prodnDR: %e prodnDL: %e\n", FunctGroupArray[guild].groupCode, cohort, initialBiomass, FunctGroupArray[guild].prodnDR[cohort], FunctGroupArray[guild].prodnDL[cohort]);
 			Group_Transfer_Contaminant(bm, boxLayerInfo, WC, WC, RefDetIndex, 0, guild, cohort, FunctGroupArray[guild].prodnDR[cohort], 0, initialBiomass, bm->dtsz_stored, 0, 0);
 			Group_Transfer_Contaminant(bm, boxLayerInfo, WC, WC, LabDetIndex, 0, guild, cohort, FunctGroupArray[guild].prodnDL[cohort] + FunctGroupArray[guild].lysis[cohort], 0, initialBiomass, bm->dtsz_stored, 0, 1);
+			/* Transfer contaminants to CARRION (dead organic matter) from total mortality */
+			/* Use total detritus production (DL + DR) as the mortality amount */
+			double totalMortality = FunctGroupArray[guild].prodnDL[cohort] + FunctGroupArray[guild].prodnDR[cohort];
+			Group_Transfer_Contaminant(bm, boxLayerInfo, WC, WC, CarrionIndex, 0, guild, cohort, totalMortality, 0, initialBiomass, bm->dtsz_stored, 0, 2);
 
 			/* Release contaminants back into the water column */
             //fprintf(bm->logFile,"Calling Gain_Contaminants from Update_Detritus WC WC - %s-%d with initialBiomass: %e releaseNH: %e\n", FunctGroupArray[guild].groupCode, cohort, initialBiomass, FunctGroupArray[guild].releaseNH[cohort]);

atlantis/atecology/atprocess.c

@@ -2262,15 +2262,11 @@ void UpdateTrackedMort(MSEBoxModel *bm, FILE *llogfp, int guildcase, int cohort,
 	}
 
 	if(bm->track_contaminants){
-
-        //fprintf(llogfp, "Time: %e box %d-%d %s-%d eaten by %s-%d habitat %d GRAZEinfo: %e scalar: %e EATINGinfo: %e\n", bm->dayt, bm->current_box, bm->current_layer, FunctGroupArray[prey].groupCode, prey_chrt, FunctGroupArray[guildcase].groupCode, cohort, habitat, GRAZEinfo[prey][prey_chrt][habitat], scalar, EATINGinfo[prey][prey_chrt][habitat]);
-
+
         if(FunctGroupArray[guildcase].isOncePerDt == TRUE) {
-
-            Group_Transfer_Contaminant(bm, boxLayerInfo, habitat, habitat, guildcase, cohort, prey, prey_chrt, (GRAZEinfo[prey][prey_chrt][habitat] * scalar), VERTinfo, EATINGinfo[prey][prey_chrt][habitat], bm->dtsz_stored, 1, 31);
+            Group_Transfer_Contaminant(bm, boxLayerInfo, level_id, habitat, guildcase, cohort, prey, prey_chrt, (GRAZEinfo[prey][prey_chrt][habitat] * scalar), VERTinfo, EATINGinfo[prey][prey_chrt][habitat], bm->dtsz_stored, 1, 31);
         } else {
-
-            Group_Transfer_Contaminant(bm, boxLayerInfo, habitat, habitat, guildcase, cohort, prey, prey_chrt, (GRAZEinfo[prey][prey_chrt][habitat] * scalar), VERTinfo, EATINGinfo[prey][prey_chrt][habitat], bm->dt, 1, 32);
+            Group_Transfer_Contaminant(bm, boxLayerInfo, level_id, habitat, guildcase, cohort, prey, prey_chrt, (GRAZEinfo[prey][prey_chrt][habitat] * scalar), VERTinfo, EATINGinfo[prey][prey_chrt][habitat], bm->dt, 1, 32);
         }
 	}
 

atlantis/atlantisUtil/atUtilGroupIO.c

@@ -664,16 +664,15 @@ int Util_Read_Functional_Group_XML(MSEBoxModel *bm, char *fileName, FILE *llogfp
 			FunctGroupArray[guild].diagTol = 2;
 		} else {
 			if (FunctGroupArray[guild].groupType == CARRION)
-				FunctGroupArray[guild].diagTol = 3;
+				FunctGroupArray[guild].diagTol = 1;  // Changed from 3 to 1 so CARRION can process multiple times per timestep for contaminant accumulation
 			else if (FunctGroupArray[guild].groupType == SED_EP_OTHER ||
 					FunctGroupArray[guild].groupType == MOB_EP_OTHER	||
-					FunctGroupArray[guild].groupType == PWN ||
 					FunctGroupArray[guild].groupType == CEP ||
 					FunctGroupArray[guild].groupType == LG_ZOO ||
 					FunctGroupArray[guild].groupType == LG_INF)
 				FunctGroupArray[guild].diagTol = 2;
 			else {
-				FunctGroupArray[guild].diagTol = 1;  // So all phytoplankton, small zooplankton and filter feeders set the adaptive time step
+				FunctGroupArray[guild].diagTol = 1;  // So all phytoplankton, small zooplankton, filter feeders, and PWN set the adaptive time step to allow multiple feeding/uptake cycles
 				if(FunctGroupArray[guild].isCultured) // Except if an aquaculture species then just have it inflicted upon them
 					FunctGroupArray[guild].diagTol = 0;
 			}

atlantis/outamps

@@ -0,0 +1 @@
+bash: cd: too many arguments