diff --git a/example-configs/skipjack/skipjack.xml b/example-configs/skipjack/skipjack.xml
index 0561453..f887e11 100755
--- a/example-configs/skipjack/skipjack.xml
+++ b/example-configs/skipjack/skipjack.xml
@@ -144,6 +144,8 @@
         <nb_files skj="22"/>
         <skj file1="skj-pop-tags_20090115.txt" file2="skj-pop-tags_20090215.txt" file3="skj-pop-tags_20090315.txt" file4="skj-pop-tags_20090415.txt" file5="skj-pop-tags_20090515.txt" file6="skj-pop-tags_20090615.txt" file7="skj-pop-tags_20090815.txt" file8="skj-pop-tags_20091015.txt" file9="skj-pop-tags_20091115.txt" file10="skj-pop-tags_20091215.txt" file11="skj-pop-tags_20100115.txt" file12="skj-pop-tags_20100215.txt" file13="skj-pop-tags_20100315.txt" file14="skj-pop-tags_20100415.txt" file15="skj-pop-tags_20100515.txt" file16="skj-pop-tags_20100615.txt" file17="skj-pop-tags_20100715.txt" file18="skj-pop-tags_20100815.txt" file19="skj-pop-tags_20100915.txt" file20="skj-pop-tags_20101015.txt" file21="skj-pop-tags_20101115.txt" file22="skj-pop-tags_20101215.txt"/>
 </file_tag_data>
+<use_tag_masks flag="0"/>
+<strdir_tagmasks value="/home/lucasb/SEAPODYM/data/tagging/skj/ds2016/tagmasks/run-interim_2x30d_po/biomass5pcts/0.5knots/"/>
 
 <!-- Output Directory and writing flags-->
 <strdir_output value="output/"/>
diff --git a/example-configs/skipjack/skipjack_F0.xml b/example-configs/skipjack/skipjack_F0.xml
index 4132140..a08e268 100755
--- a/example-configs/skipjack/skipjack_F0.xml
+++ b/example-configs/skipjack/skipjack_F0.xml
@@ -144,6 +144,8 @@
         <nb_files skj="22"/>
         <skj file1="skj-pop-tags_20090115.txt" file2="skj-pop-tags_20090215.txt" file3="skj-pop-tags_20090315.txt" file4="skj-pop-tags_20090415.txt" file5="skj-pop-tags_20090515.txt" file6="skj-pop-tags_20090615.txt" file7="skj-pop-tags_20090815.txt" file8="skj-pop-tags_20091015.txt" file9="skj-pop-tags_20091115.txt" file10="skj-pop-tags_20091215.txt" file11="skj-pop-tags_20100115.txt" file12="skj-pop-tags_20100215.txt" file13="skj-pop-tags_20100315.txt" file14="skj-pop-tags_20100415.txt" file15="skj-pop-tags_20100515.txt" file16="skj-pop-tags_20100615.txt" file17="skj-pop-tags_20100715.txt" file18="skj-pop-tags_20100815.txt" file19="skj-pop-tags_20100915.txt" file20="skj-pop-tags_20101015.txt" file21="skj-pop-tags_20101115.txt" file22="skj-pop-tags_20101215.txt"/>
 </file_tag_data>
+<use_tag_masks flag="0"/>
+<strdir_tagmasks value="/home/lucasb/SEAPODYM/data/tagging/skj/ds2016/tagmasks/run-interim_2x30d_po/biomass5pcts/0.5knots/"/>
 
 <!-- Output Directory and writing flags-->
 <strdir_output value="output/"/>
diff --git a/src/Map.cpp b/src/Map.cpp
index 948dc0b..ebcd520 100755
--- a/src/Map.cpp
+++ b/src/Map.cpp
@@ -1,126 +1,187 @@
 #include "Map.h"
 
 ////////////////////////////////////////////////////////////
-// fonction lit_map de la classe PMap 
+// fonction lit_map de la classe PMap
 // lit le fichier source et charge les donnees dans les tableaux
 
 int nbl;
 double deltaX;
-//CBord bord_layer;
-
+// CBord bord_layer;
 
 void PMap::lit_map(CParam &param)
 {
-	const int nti=param.get_nbi();
-	const int ntj=param.get_nbj();
-	nbl=param.nb_layer;
-
-// creation des tableaux
+	const int nti = param.get_nbi();
+	const int ntj = param.get_nbj();
+	nbl = param.nb_layer;
 
+	// creation des tableaux
 
 	bord_cell.allocate(0, nti - 1, 0, ntj - 1);
 	nbl_bord_cell.allocate(0, nti - 1, 0, ntj - 1);
-	carte.allocate(0, nti - 1, 0, ntj - 1);// carte
-	itopo.allocate(0, nti - 1, 0, ntj - 1);// carte
+	carte.allocate(0, nti - 1, 0, ntj - 1); // carte
+	itopo.allocate(0, nti - 1, 0, ntj - 1); // carte
 	carte.initialize();
-	itopo.initialize(); itopo = 1; //default value
+	itopo.initialize();
+	itopo = 1; // default value
 
-	//indices for regions
+	// indices for regions
 	reg_indices(param);
-	
-// lecture du fichier
-	
+
+	// lecture du fichier
+
 	ifstream litcarte(param.str_file_mask.c_str(), ios::in);
 
-	if (!litcarte) {cout<<endl<<"WARNING : Cannot read file: "<< param.str_file_mask.c_str() << endl;}
-	for (int j=1; j<ntj-1; j++)
+	if (!litcarte)
 	{
-		for (int i=1; i<nti-1; i++)
+		cout << endl
+			 << "WARNING : Cannot read file: " << param.str_file_mask.c_str() << endl;
+	}
+	for (int j = 1; j < ntj - 1; j++)
+	{
+		for (int i = 1; i < nti - 1; i++)
 		{
 			litcarte >> carte[i][j];
-			if (carte[i][j]>nbl) carte[i][j] = nbl;
-
+			if (carte[i][j] > nbl)
+				carte[i][j] = nbl;
 		}
 	}
 	litcarte.close();
 
 	ifstream littopo(param.str_file_topo.c_str(), ios::in);
 
-	if (!littopo) {cout<<endl<<"WARNING : Cannot read file: "<< param.str_file_topo.c_str() << endl;}
-	for (int j=1; j<ntj-1; j++)
+	if (!littopo)
 	{
-		for (int i=1; i<nti-1; i++)
+		cout << endl
+			 << "WARNING : Cannot read file: " << param.str_file_topo.c_str() << endl;
+	}
+	for (int j = 1; j < ntj - 1; j++)
+	{
+		for (int i = 1; i < nti - 1; i++)
 		{
 			littopo >> itopo[i][j];
-			//necessary to avoid Ha=0
-			if (carte[i][j] && itopo(i,j)==0) itopo(i,j)=0.0001;
+			// necessary to avoid Ha=0
+			if (carte[i][j] && itopo(i, j) == 0)
+				itopo(i, j) = 0.0001;
 		}
 	}
 	littopo.close();
 
-	if(param.nb_EEZ)
+	if (param.nb_EEZ)
 	{
 		// mask contenant les codes des Zones exclusives
-		maskEEZ.allocate(0, nti-1, 0, ntj-1);
+		maskEEZ.allocate(0, nti - 1, 0, ntj - 1);
 		maskEEZ.initialize();
 
 		ifstream rtxt(param.str_file_maskEEZ.c_str(), ios::in);
-		if (!rtxt) {cout<<endl<<"WARNING : Cannot read file: "<< 
-			param.str_file_maskEEZ.c_str() << endl;}
+		if (!rtxt)
+		{
+			cout << endl
+				 << "WARNING : Cannot read file: " << param.str_file_maskEEZ.c_str() << endl;
+		}
 
-		for (int j=1;j<ntj-1;j++)
-			for (int i=1;i<nti-1;i++)
+		for (int j = 1; j < ntj - 1; j++)
+			for (int i = 1; i < nti - 1; i++)
 			{
 				rtxt >> maskEEZ[i][j];
 			}
-		rtxt.close(); 
+		rtxt.close();
 	}
 
-	if(param.mpa_simulation)
+	if (param.mpa_simulation)
 	{
 		// mask contenant les codes des Zones exclusives
-		maskMPA.allocate(0, nti-1, 0, ntj-1);
+		maskMPA.allocate(0, nti - 1, 0, ntj - 1);
 		maskMPA.initialize();
 
 		ifstream rtxt(param.str_file_maskMPA.c_str(), ios::in);
-		if (!rtxt) {cout<<endl<<"WARNING : Cannot read file: "<< param.str_file_maskMPA.c_str() << endl;}
-		for (int j=1;j<ntj-1;j++)
-			for (int i=1;i<nti-1;i++)
+		if (!rtxt)
+		{
+			cout << endl
+				 << "WARNING : Cannot read file: " << param.str_file_maskMPA.c_str() << endl;
+		}
+		for (int j = 1; j < ntj - 1; j++)
+			for (int i = 1; i < nti - 1; i++)
 			{
 				rtxt >> maskMPA[i][j];
 			}
-		rtxt.close(); 
+		rtxt.close();
 	}
 
-
-	//initialize
-	global = 0;//Attn: for optimization no global until fixed!
+	// initialize
+	global = 0; // Attn: for optimization no global until fixed!
 	deltaX = param.deltaX;
-	double reso_x = param.deltaX/60.0;
-	domain_type((nti-2)*reso_x);
-//global = 0;//ATTN: put global to zero to be able to run with summy fisheries. Re-set when needed.
-	definit_cell_bords(nti,ntj); // appel de la fonction 
-	definit_lim_infsup(nti,ntj); // appel de la fonction
+	double reso_x = param.deltaX / 60.0;
+	domain_type((nti - 2) * reso_x);
+	// global = 0;//ATTN: put global to zero to be able to run with summy fisheries. Re-set when needed.
+	definit_cell_bords(nti, ntj); // appel de la fonction
+	definit_lim_infsup(nti, ntj); // appel de la fonction
 }
 
-void PMap::reg_indices(CParam& param)
+void PMap::lit_tagmap(CParam &param, int tagpop)
 {
-	//Regional indices
-	const int nb_regs = param.nb_region;
-	regimin.allocate(0,nb_regs-1);
-	regjmin.allocate(0,nb_regs-1);
-	regimax.allocate(0,nb_regs-1);
-	regjmax.allocate(0,nb_regs-1);
-	for (int a =0; a < param.nb_region; a++){
+	const int nti = param.get_nbi();
+	const int ntj = param.get_nbj();
+	nbl = param.nb_layer;
 
-		double lonmin = param.area[a]->lgmin; double lonmax = param.area[a]->lgmax;
-		double latmin = param.area[a]->ltmax; double latmax = param.area[a]->ltmin;
-		regimin[a] = param.lontoi(lonmin); regimax[a] = param.lontoi(lonmax);
-		regjmin[a] = param.lattoj(latmin); regjmax[a] = param.lattoj(latmax);
+	// creation des tableaux
 
+	bord_cell.allocate(0, nti - 1, 0, ntj - 1);
+	nbl_bord_cell.allocate(0, nti - 1, 0, ntj - 1);
+	carte.allocate(0, nti - 1, 0, ntj - 1); // carte
+	carte.initialize();
 
-//cout << a+1 << "\t" << param.area[a]->lgmin << "\t" << param.area[a]->lgmax << "\t" << param.area[a]->ltmin << "\t" << param.area[a]->ltmax << endl;
-//cout << a+1 << "\t" << regimin[a] << "\t" << regimax[a] << "\t" << regjmin[a] << "\t" << regjmax[a] << endl;
+	// lecture du fichier
+
+	ifstream litcarte(param.file_tag_masks[tagpop].c_str(), ios::in);
+
+	if (!litcarte)
+	{
+		cout << endl
+			 << "WARNING : Cannot read file: " << param.str_file_mask.c_str() << endl;
+	}
+	for (int j = 1; j < ntj - 1; j++)
+	{
+		for (int i = 1; i < nti - 1; i++)
+		{
+			litcarte >> carte[i][j];
+			if (carte[i][j] > nbl)
+				carte[i][j] = nbl;
+		}
+	}
+	litcarte.close();
+
+	// initialize
+	global = 0; // Attn: for optimization no global until fixed!
+	deltaX = param.deltaX;
+	double reso_x = param.deltaX / 60.0;
+	domain_type((nti - 2) * reso_x);
+	// global = 0;//ATTN: put global to zero to be able to run with summy fisheries. Re-set when needed.
+	definit_cell_bords(nti, ntj); // appel de la fonction
+	definit_lim_infsup(nti, ntj); // appel de la fonction
+}
+
+void PMap::reg_indices(CParam &param)
+{
+	// Regional indices
+	const int nb_regs = param.nb_region;
+	regimin.allocate(0, nb_regs - 1);
+	regjmin.allocate(0, nb_regs - 1);
+	regimax.allocate(0, nb_regs - 1);
+	regjmax.allocate(0, nb_regs - 1);
+	for (int a = 0; a < param.nb_region; a++)
+	{
+
+		double lonmin = param.area[a]->lgmin;
+		double lonmax = param.area[a]->lgmax;
+		double latmin = param.area[a]->ltmax;
+		double latmax = param.area[a]->ltmin;
+		regimin[a] = param.lontoi(lonmin);
+		regimax[a] = param.lontoi(lonmax);
+		regjmin[a] = param.lattoj(latmin);
+		regjmax[a] = param.lattoj(latmax);
+
+		// cout << a+1 << "\t" << param.area[a]->lgmin << "\t" << param.area[a]->lgmax << "\t" << param.area[a]->ltmin << "\t" << param.area[a]->ltmax << endl;
+		// cout << a+1 << "\t" << regimin[a] << "\t" << regimax[a] << "\t" << regjmin[a] << "\t" << regjmax[a] << endl;
 	}
 }
 /*
@@ -128,587 +189,664 @@ char PMap::get_bord_layer_x(const int i, const int j){
 
 	bord_layer.b	= nbl_bord_cell[i][j];
 	char pos	= bord_layer.cotex();
-	return pos;	
+	return pos;
 }
 
 char PMap::get_bord_layer_y(const int i, const int j){
 
 	bord_layer.b	= nbl_bord_cell[i][j];
 	char pos	= bord_layer.cotey();
-	return pos;	
+	return pos;
 }
 */
-void PMap::definit_cell_bords(const int nti, const int ntj) // définition des bords des cellules de la zone
+void PMap::definit_cell_bords(const int nti, const int ntj) // définition des bords des cellules de la zone
 {
-  CBord bordures; //bordures, objet de la classe Cbords  
-
-  int i, j;
-  
-  // Initialisation des cellules avec tous les cotes ouverts
-  // SANS = sans frontiere, tous cotes ouverts
-  // G_FERME = ferme a gauche
-  // D_FERME = ferme a droite
-
-  bordures.b = SANS;
-  for (i = 1; i < nti-1; i++)
-  {
-    for (j = 1; j < ntj-1; j++)
-    {
-      bord_cell[i][j] = bordures.b;
-      nbl_bord_cell[i][j] = bordures.b;
-    }
-  }  
-  //***********************
-  // Cadre de la grille 
-  // Definition des cotés des cellules situées aux bords de la grille
-  //***********************
-  // bord horizontal haut de la grille, donc ferme a gauche
-  i = 1;
-  bordures.cote.y = SANS;
-  bordures.cote.x = G_FERME;
-//cout << "i=1 " << (int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl;
-  for (j = 1; j < ntj-1; j++)
-  {
-    bord_cell[i][j] = bordures.b;
-    nbl_bord_cell[i][j] = bordures.b;
-  }
-  
-  // bord horizontal bas de la grille donc ferme a droite
-  i = nti-2;
-  bordures.cote.y = SANS;
-  bordures.cote.x = D_FERME;
-//cout << "i=nbi-2 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " " << bordures.b << endl;
-  for (j = 1; j < ntj-1; j++)
-  {
-    bord_cell[i][j] = bordures.b;
-    nbl_bord_cell[i][j] = bordures.b;
-  }
-
-  // bord vertical gauche de la grille donc ferme a gauche
-  j = 1;
-  bordures.cote.x = SANS;
-  bordures.cote.y = G_FERME;
-//cout << "j=1 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl;
-  for (i = 1; i < nti-1; i++)
-  {
-    bord_cell[i][j] = bordures.b;
-    nbl_bord_cell[i][j] = bordures.b;
-  }
-
-  // bord vertical droit de la grille donc ferme a droite
-  j = ntj-2;
-  bordures.cote.x = SANS;
-  bordures.cote.y = D_FERME;
-//cout << "j=nbj-2 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl;
-  for (i = 1; i < nti-1; i++)
-  {
-    bord_cell[i][j] = bordures.b;
-    nbl_bord_cell[i][j] = bordures.b;
-  }
-
-  // finir les coins du cadre de la grille
-  
-  //coin en haut a gauche  
-  bordures.cote.x = G_FERME;
-  bordures.cote.y = G_FERME;
-  bord_cell[1][1] = bordures.b;
-  nbl_bord_cell[1][1] = bordures.b;
-//cout << "i=1; j=1 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl;
-  
-  //coin en haut a droite
-  bordures.cote.x = G_FERME;
-  bordures.cote.y = D_FERME;
-  bord_cell[1][ntj-2]=bordures.b;
-  nbl_bord_cell[1][ntj-2]=bordures.b;
-//cout << "i=1; j=nbj-2 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl;
-  
-  // coin en bas a gauche
-  bordures.cote.x = D_FERME;
-  bordures.cote.y = G_FERME;
-  bord_cell[nti-2][1] = bordures.b;
-  nbl_bord_cell[nti-2][1] = bordures.b;
-//cout << "i=nbi-2; j=1 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl;
-
-  // coin en bas a droite
-  bordures.cote.x = D_FERME;
-  bordures.cote.y = D_FERME;
-  bord_cell[nti-2][ntj-2] = bordures.b; 
-  nbl_bord_cell[nti-2][ntj-2] = bordures.b; 
-//cout << "i=nbi-2; j=nbj-2 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl; 
-//exit(1);
-  //************************
-  // Interieur de la grille
-  // cas ou il y a des terres isolees 
-  //************************
-  for (i = 1; i < nti-1; i++)  {
-    for (j = 1; j < ntj-1; j++){
-      if (!carte[i][j]) {// si la valeur de la carte = 0 (donc terre)
-		if (i>1){	// si cellule n'est pas sur le bord horizontal haut
-		
-			bordures.b = bord_cell[i-1][j];
-			if (bordures.cotex() == SANS)
-			{
-			bordures.cote.x = D_FERME;
-			bord_cell[i-1][j] = bordures.b;
-			}
+	CBord bordures; // bordures, objet de la classe Cbords
+
+	int i, j;
+
+	// Initialisation des cellules avec tous les cotes ouverts
+	// SANS = sans frontiere, tous cotes ouverts
+	// G_FERME = ferme a gauche
+	// D_FERME = ferme a droite
+
+	bordures.b = SANS;
+	for (i = 1; i < nti - 1; i++)
+	{
+		for (j = 1; j < ntj - 1; j++)
+		{
+			bord_cell[i][j] = bordures.b;
+			nbl_bord_cell[i][j] = bordures.b;
 		}
-		if (i<nti-2){ // si cellule n'est pas sur le bord horizontal bas
-		
-			bordures.b = bord_cell[i+1][j];
-			if (bordures.cotex() == SANS)
-			{
-			bordures.cote.x = G_FERME;
-			bord_cell[i+1][j] = bordures.b;
+	}
+	//***********************
+	// Cadre de la grille
+	// Definition des cotés des cellules situées aux bords de la grille
+	//***********************
+	// bord horizontal haut de la grille, donc ferme a gauche
+	i = 1;
+	bordures.cote.y = SANS;
+	bordures.cote.x = G_FERME;
+	// cout << "i=1 " << (int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl;
+	for (j = 1; j < ntj - 1; j++)
+	{
+		bord_cell[i][j] = bordures.b;
+		nbl_bord_cell[i][j] = bordures.b;
+	}
+
+	// bord horizontal bas de la grille donc ferme a droite
+	i = nti - 2;
+	bordures.cote.y = SANS;
+	bordures.cote.x = D_FERME;
+	// cout << "i=nbi-2 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " " << bordures.b << endl;
+	for (j = 1; j < ntj - 1; j++)
+	{
+		bord_cell[i][j] = bordures.b;
+		nbl_bord_cell[i][j] = bordures.b;
+	}
+
+	// bord vertical gauche de la grille donc ferme a gauche
+	j = 1;
+	bordures.cote.x = SANS;
+	bordures.cote.y = G_FERME;
+	// cout << "j=1 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl;
+	for (i = 1; i < nti - 1; i++)
+	{
+		bord_cell[i][j] = bordures.b;
+		nbl_bord_cell[i][j] = bordures.b;
+	}
+
+	// bord vertical droit de la grille donc ferme a droite
+	j = ntj - 2;
+	bordures.cote.x = SANS;
+	bordures.cote.y = D_FERME;
+	// cout << "j=nbj-2 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl;
+	for (i = 1; i < nti - 1; i++)
+	{
+		bord_cell[i][j] = bordures.b;
+		nbl_bord_cell[i][j] = bordures.b;
+	}
+
+	// finir les coins du cadre de la grille
+
+	// coin en haut a gauche
+	bordures.cote.x = G_FERME;
+	bordures.cote.y = G_FERME;
+	bord_cell[1][1] = bordures.b;
+	nbl_bord_cell[1][1] = bordures.b;
+	// cout << "i=1; j=1 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl;
+
+	// coin en haut a droite
+	bordures.cote.x = G_FERME;
+	bordures.cote.y = D_FERME;
+	bord_cell[1][ntj - 2] = bordures.b;
+	nbl_bord_cell[1][ntj - 2] = bordures.b;
+	// cout << "i=1; j=nbj-2 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl;
+
+	// coin en bas a gauche
+	bordures.cote.x = D_FERME;
+	bordures.cote.y = G_FERME;
+	bord_cell[nti - 2][1] = bordures.b;
+	nbl_bord_cell[nti - 2][1] = bordures.b;
+	// cout << "i=nbi-2; j=1 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl;
+
+	// coin en bas a droite
+	bordures.cote.x = D_FERME;
+	bordures.cote.y = D_FERME;
+	bord_cell[nti - 2][ntj - 2] = bordures.b;
+	nbl_bord_cell[nti - 2][ntj - 2] = bordures.b;
+	// cout << "i=nbi-2; j=nbj-2 " <<(int)bordures.cote.x <<" " <<(int)bordures.cote.y << " "<< bordures.b << endl;
+	// exit(1);
+	//************************
+	//  Interieur de la grille
+	//  cas ou il y a des terres isolees
+	//************************
+	for (i = 1; i < nti - 1; i++)
+	{
+		for (j = 1; j < ntj - 1; j++)
+		{
+			if (!carte[i][j])
+			{ // si la valeur de la carte = 0 (donc terre)
+				if (i > 1)
+				{ // si cellule n'est pas sur le bord horizontal haut
+
+					bordures.b = bord_cell[i - 1][j];
+					if (bordures.cotex() == SANS)
+					{
+						bordures.cote.x = D_FERME;
+						bord_cell[i - 1][j] = bordures.b;
+					}
+				}
+				if (i < nti - 2)
+				{ // si cellule n'est pas sur le bord horizontal bas
+
+					bordures.b = bord_cell[i + 1][j];
+					if (bordures.cotex() == SANS)
+					{
+						bordures.cote.x = G_FERME;
+						bord_cell[i + 1][j] = bordures.b;
+					}
+				}
+				if (j > 1)
+				{ // si cellule n'est pas sur le bord vertical gauche
+
+					bordures.b = bord_cell[i][j - 1];
+					if (bordures.cotey() == SANS)
+					{
+						bordures.cote.y = D_FERME;
+						bord_cell[i][j - 1] = bordures.b;
+					}
+				}
+				if (j < ntj - 2)
+				{ // si cellule n'est pas sur le bord vertical droit
+
+					bordures.b = bord_cell[i][j + 1];
+					if (bordures.cotey() == SANS)
+					{
+						bordures.cote.y = G_FERME;
+						bord_cell[i][j + 1] = bordures.b;
+					}
+				}
+				bordures.cote.x = TERRE;
+				bordures.cote.y = TERRE;
+				bord_cell[i][j] = bordures.b;
 			}
-		}
-		if (j>1){ // si cellule n'est pas sur le bord vertical gauche
-		
-			bordures.b = bord_cell[i][j-1];
-			if (bordures.cotey() == SANS)
-			{
-			bordures.cote.y = D_FERME;
-			bord_cell[i][j-1] = bordures.b;
+		} // fin de boucle j
+	}	  // fin de boucle i
+
+	// set the boundaries for nbl-layer cells
+	for (i = 1; i < nti - 1; i++)
+	{
+		for (j = 1; j < ntj - 1; j++)
+		{
+			if (carte[i][j] < nbl)
+			{ // if number of layers in the cell < total nbl
+				if (i > 1)
+				{ // si cellule n'est pas sur le bord horizontal haut
+
+					bordures.b = nbl_bord_cell[i - 1][j];
+					if (bordures.cotex() == SANS)
+					{
+						bordures.cote.x = D_FERME;
+						nbl_bord_cell[i - 1][j] = bordures.b;
+					}
+				}
+				if (i < nti - 2)
+				{ // si cellule n'est pas sur le bord horizontal bas
+
+					bordures.b = nbl_bord_cell[i + 1][j];
+					if (bordures.cotex() == SANS)
+					{
+						bordures.cote.x = G_FERME;
+						nbl_bord_cell[i + 1][j] = bordures.b;
+					}
+				}
+				if (j > 1)
+				{ // si cellule n'est pas sur le bord vertical gauche
+
+					bordures.b = nbl_bord_cell[i][j - 1];
+					if (bordures.cotey() == SANS)
+					{
+						bordures.cote.y = D_FERME;
+						nbl_bord_cell[i][j - 1] = bordures.b;
+					}
+				}
+				if (j < ntj - 2)
+				{ // si cellule n'est pas sur le bord vertical droit
+
+					bordures.b = nbl_bord_cell[i][j + 1];
+					if (bordures.cotey() == SANS)
+					{
+						bordures.cote.y = G_FERME;
+						nbl_bord_cell[i][j + 1] = bordures.b;
+					}
+				}
+				bordures.cote.x = TERRE;
+				bordures.cote.y = TERRE;
+				nbl_bord_cell[i][j] = bordures.b;
 			}
+		} // fin de boucle j
+	}	  // fin de boucle i
+
+	//******************
+	// cas ou il ya des Terres dans les bords de la grille
+	//******************
+	i = 1;
+	for (j = 1; j < ntj - 1; j++)
+	{
+		bordures.b = bord_cell[i][j];
+		if (carte[i + 1][j] == 0)
+		{
+			bordures.cote.x = TERRE;
+			bord_cell[i][j] = bordures.b;
 		}
-		if (j<ntj-2){ // si cellule n'est pas sur le bord vertical droit
-		
-			bordures.b = bord_cell[i][j+1];
-			if (bordures.cotey() == SANS)
-			{
-		bordures.cote.y = G_FERME;
-		bord_cell[i][j+1] = bordures.b;
-		  	}
+	}
+
+	i = nti - 2;
+	for (j = 1; j < ntj - 1; j++)
+	{
+		bordures.b = bord_cell[i][j];
+		if (carte[i - 1][j] == 0)
+		{
+			bordures.cote.x = TERRE;
+			bord_cell[i][j] = bordures.b;
 		}
-		bordures.cote.x = TERRE;
-		bordures.cote.y = TERRE;
-		bord_cell[i][j] = bordures.b;
-      }
-	} // fin de boucle j
-  } // fin de boucle i
-
-  //set the boundaries for nbl-layer cells
-  for (i = 1; i < nti-1; i++)  {
-    for (j = 1; j < ntj-1; j++){
-      if (carte[i][j]<nbl) { // if number of layers in the cell < total nbl
-		if (i>1){	// si cellule n'est pas sur le bord horizontal haut
-		
-			  bordures.b = nbl_bord_cell[i-1][j];
-			  if (bordures.cotex() == SANS){
-			    bordures.cote.x = D_FERME;
-			    nbl_bord_cell[i-1][j] = bordures.b;
-			  }
+	}
+
+	j = 1;
+	for (i = 1; i < nti - 1; i++)
+	{
+		bordures.b = bord_cell[i][j];
+		if (carte[i][j + 1] == 0)
+		{
+			bordures.cote.y = TERRE;
+			bord_cell[i][j] = bordures.b;
 		}
-		if (i<nti-2){ // si cellule n'est pas sur le bord horizontal bas
-		
-			bordures.b = nbl_bord_cell[i+1][j];
-			if (bordures.cotex() == SANS){
-			  bordures.cote.x = G_FERME;
-			  nbl_bord_cell[i+1][j] = bordures.b;
-			}
+	}
+
+	j = ntj - 2;
+	for (i = 1; i < nti - 1; i++)
+	{
+		bordures.b = bord_cell[i][j];
+		if (carte[i][j - 1] == 0)
+		{
+			bordures.cote.y = TERRE;
+			bord_cell[i][j] = bordures.b;
 		}
-		if (j>1){ // si cellule n'est pas sur le bord vertical gauche
-		
-			bordures.b = nbl_bord_cell[i][j-1];
-			if (bordures.cotey() == SANS){
-			  bordures.cote.y = D_FERME;
-			  nbl_bord_cell[i][j-1] = bordures.b;
-			}
+	}
+
+	// once more, similarly let's set boundaries for nb_layer cells
+	i = 1;
+	for (j = 1; j < ntj - 1; j++)
+	{
+		bordures.b = nbl_bord_cell[i][j];
+		if (carte[i + 1][j] < nbl)
+		{
+			bordures.cote.x = TERRE;
+			nbl_bord_cell[i][j] = bordures.b;
 		}
-		if (j<ntj-2){ // si cellule n'est pas sur le bord vertical droit
-		
-			bordures.b = nbl_bord_cell[i][j+1];
-			if (bordures.cotey() == SANS){
-	   		  bordures.cote.y = G_FERME;
-			  nbl_bord_cell[i][j+1] = bordures.b;
-		  	}
+	}
+
+	i = nti - 2;
+	for (j = 1; j < ntj - 1; j++)
+	{
+		bordures.b = nbl_bord_cell[i][j];
+		if (carte[i - 1][j] < nbl)
+		{
+			bordures.cote.x = TERRE;
+			nbl_bord_cell[i][j] = bordures.b;
 		}
-		bordures.cote.x = TERRE;
-		bordures.cote.y = TERRE;
-		nbl_bord_cell[i][j] = bordures.b;
-      }
-	} // fin de boucle j
-  } // fin de boucle i
-
-
-  //******************
-  // cas ou il ya des Terres dans les bords de la grille
-  //******************
-  i=1;
-  		for (j=1 ; j < ntj-1 ; j++){
-  			bordures.b=bord_cell[i][j];
-  			if (carte[i+1][j]==0) {
-	  			bordures.cote.x=TERRE;
-	  			bord_cell[i][j]=bordures.b;
+	}
+	// Inna Oct05: open left-right boundary
+	if (global)
+	{
+		i = 1;
+		for (j = 1; j < ntj - 1; j++)
+		{
+			bordures.b = bord_cell[i][j];
+			if (carte[1][j] && carte[nti - 2][j])
+			{
+				bordures.cote.x = SANS;
+				bord_cell[i][j] = bordures.b;
 			}
 		}
-		
-  i=nti-2;
-  		for (j=1 ; j < ntj-1 ; j++){
-  			bordures.b=bord_cell[i][j];
-  			if (carte[i-1][j]==0) {
-	  			bordures.cote.x=TERRE;
-	  			bord_cell[i][j]=bordures.b;
-	  		}
-  		}
-
-  j=1;
-  		for (i=1 ; i < nti-1 ; i++){
-  			bordures.b=bord_cell[i][j];
-  			if (carte[i][j+1]==0) {
-	  			bordures.cote.y=TERRE;
-	  			bord_cell[i][j]=bordures.b;
-	  		}
-	  	}
-	  	
-  j=ntj-2;
-  		for (i=1 ; i < nti-1 ; i++){
-  			bordures.b=bord_cell[i][j];
-  			if (carte[i][j-1]==0) {
-	  			bordures.cote.y=TERRE;
-	  			bord_cell[i][j]=bordures.b;
-	  		}
-	  	}
-
-  //once more, similarly let's set boundaries for nb_layer cells
-  i=1;
-  		for (j=1 ; j < ntj-1 ; j++){
-  			bordures.b=nbl_bord_cell[i][j];
-  			if (carte[i+1][j]<nbl) {
-	  			bordures.cote.x=TERRE;
-	  			nbl_bord_cell[i][j]=bordures.b;
+
+		i = nti - 2;
+		for (j = 1; j < ntj - 1; j++)
+		{
+			bordures.b = bord_cell[i][j];
+			if (carte[1][j] && carte[nti - 2][j])
+			{
+				bordures.cote.x = SANS;
+				bord_cell[i][j] = bordures.b;
 			}
 		}
-		
-  i=nti-2;
-  		for (j=1 ; j < ntj-1 ; j++){
-  			bordures.b=nbl_bord_cell[i][j];
-  			if (carte[i-1][j]<nbl) {
-	  			bordures.cote.x=TERRE;
-	  			nbl_bord_cell[i][j]=bordures.b;
-	  		}
-  		}
-        //Inna Oct05: open left-right boundary
-        if (global){
-                i = 1;
-                for (j=1 ; j < ntj-1 ; j++){
-                        bordures.b=bord_cell[i][j];
-                        if (carte[1][j] && carte[nti-2][j]){
-                                bordures.cote.x=SANS;
-                                bord_cell[i][j]=bordures.b;
-                        }
-                }
-
-                i=nti-2;
-                for (j=1 ; j < ntj-1 ; j++){
-                        bordures.b=bord_cell[i][j];
-                        if (carte[1][j] && carte[nti-2][j])
-                        {
-                                bordures.cote.x=SANS;
-                                bord_cell[i][j]=bordures.b;
-                        }
-                }
-        }
-//END OF TEST		
-		
-
-  j=1;
-  		for (i=1 ; i < nti-1 ; i++){
-  			bordures.b=nbl_bord_cell[i][j];
-  			if (carte[i][j+1]<nbl) {
-	  			bordures.cote.y=TERRE;
-	  			nbl_bord_cell[i][j]=bordures.b;
-	  		}
-	  	}
-
-	  	
-  j=ntj-2;
-  		for (i=1 ; i < nti-1 ; i++){
-  			bordures.b=nbl_bord_cell[i][j];
-  			if (carte[i][j-1]<nbl) {
-	  			bordures.cote.y=TERRE;
-	  			nbl_bord_cell[i][j]=bordures.b;
-	  		}
-	  	}
-}   
+	}
+	// END OF TEST
+
+	j = 1;
+	for (i = 1; i < nti - 1; i++)
+	{
+		bordures.b = nbl_bord_cell[i][j];
+		if (carte[i][j + 1] < nbl)
+		{
+			bordures.cote.y = TERRE;
+			nbl_bord_cell[i][j] = bordures.b;
+		}
+	}
 
+	j = ntj - 2;
+	for (i = 1; i < nti - 1; i++)
+	{
+		bordures.b = nbl_bord_cell[i][j];
+		if (carte[i][j - 1] < nbl)
+		{
+			bordures.cote.y = TERRE;
+			nbl_bord_cell[i][j] = bordures.b;
+		}
+	}
+}
 
 //***********************************************
-// determine les valeurs d'indices a partir desquels 
+// determine les valeurs d'indices a partir desquels
 // il n'y a plus de terre en partant des bords de la grille
-// But: Optimization du calcul car inutile de calculer 
+// But: Optimization du calcul car inutile de calculer
 // le transport dans des cellules terre
 //************************************************
 void PMap::definit_lim_infsup(const int nti, const int ntj)
 
 {
-  const int rowmin = carte.rowmin();
-  const int rowmax = carte.rowmax();
-  const int colmin = carte.colmin();
-  const int colmax = carte.colmax();
-
-//cout << "rowmin " << rowmin << endl;
-//cout << "rowmax " << rowmax << endl;
-//cout << "colmin " << colmin << endl;
-//cout << "colmax " << colmax << endl;
-
-  IVECTOR iinftmp(colmin, colmax);  
-  IVECTOR isuptmp(colmin, colmax);  
-  iinftmp.initialize();
-  isuptmp.initialize();
-
-//cout << "iinftmp " << iinftmp << endl;
-//cout << "isuptmp " << isuptmp << endl;
-
-  for (int j = colmin; j <= colmax; j++)
-  {
-    int i = rowmin;
-    while (i <= rowmax && carte[i][j] == 0)
-      { i++; }
-    iinftmp[j] = i;
-  }
-
-//cout << "iinftmp " << iinftmp << endl;
-
-  for (int j = colmax; j >= colmin; j--)
-  {
-    int i = rowmax;
-    while (i >= rowmin && carte[i][j] == 0)
-      { i--; }
-    isuptmp[j] = i;
-  }
-
-//cout << "isuptmp" << isuptmp << endl;
-
-  jmin = colmin;
-  while (iinftmp(jmin) > isuptmp(jmin)) jmin++;
-
-  jmax = colmax;
-  while (iinftmp(jmax) > isuptmp(jmax)) jmax--;
-
-//cout << "jmin,jmax " << jmin << ' ' << jmax << endl;
-
- for (int j = jmin; j <= jmax; j++)
-  {
-	int i = iinftmp[j];
-	if (!(carte[i - 1][j] == 0 && carte[i][j] != 0))
-	{
-		//cout << carte[i][j] << ' ' << carte[i + 1][j] << endl;
-		//cout << i << ' ' << j << endl;
-		//cout << __FILE__ << ':' << __LINE__ << endl;
-		exit(1);
-	}
-
-	i = isuptmp[j];
-	if (!(carte[i][j] != 0 && carte[i + 1][j] == 0))
-	{
-		//cout << carte[i][j] << ' ' << carte[i + 1][j] << endl;
-		//cout << i << ' ' << j << endl;
-		//cout << __FILE__ << ':' << __LINE__ << endl;
-		exit(1);
-	}
-
-  }
-
-  iinf.allocate(jmin, jmax);
-  isup.allocate(jmin, jmax);
-  iinf.initialize();
-  isup.initialize();
-
-  for (int j = jmin; j <= jmax; j++)
-  {
-	  iinf[j] = iinftmp[j];
-	  isup[j] = isuptmp[j];
-  }
-//cout << "iinf " << iinf << endl;
-//cout << "isup " << isup << endl;
-//exit(1);
-//=====================================
-
-  IVECTOR jinftmp(rowmin, rowmax);  
-  IVECTOR jsuptmp(rowmin, rowmax);  
-  jinftmp.initialize();
-  jsuptmp.initialize();
-//cout << "jinftmp " << jinftmp << endl;
-//cout << "jsuptmp " << jsuptmp << endl;
-
-  for (int i = rowmin; i <= rowmax; i++)
-  {
-    int j = colmin;
-    while (j <= colmax && carte[i][j] == 0)
-      { j++; }
-    jinftmp[i] = j;
-  }
-//cout << "jinftmp " << jinftmp << endl;
-//exit(1);
-
-  for (int i = rowmax; i >= rowmin; i--)
-  {
-    int j = colmax;
-    while (j >= colmin && carte[i][j] == 0)
-      { j--; }
-    jsuptmp[i] = j;
-  }
-//cout << "jsuptmp " << jsuptmp << endl;
-//exit(1);
-
-  imin = rowmin;
-  while (jinftmp(imin) > jsuptmp(imin)) imin++;
-  imax = rowmax;
-  while (jinftmp(imax) > jsuptmp(imax)) imax--;
-//cout << imin << endl;
-//cout << imax << endl;
-//exit(1);
-
-  jinf.allocate(imin, imax);
-  jsup.allocate(imin, imax);
-  for (int i = imin; i <= imax; i++)
-  {
-    jinf[i] = jinftmp[i];
-    jsup[i] = jsuptmp[i];
-  }
-//cout << "jinf " << jinf << endl;
-//cout << "jsup " << jsup << endl;
-
-  for (int j=jmin; j <= jmax; j++)
-  {
-    int lb = iinf[j];
-    int ub = isup[j];
-    for (int i = lb; i <= ub; i++)
-    {
-      if (jinf[i] > j)
-      {
-        jinf[i] = j;
-      }
-      if (jsup[i] < j)
-      {
-	jsup[i] = j; 
-      }
-    }
-  }
-
-//cout << "jinf " << jinf << endl;
-//cout << "jsup" << jsup << endl;
-//exit(1);
-
-  for (int i = imin; i <= imax; i++)
-  {
-    int lb = jinf[i];
-    int ub = jsup[i];
-    for (int j = lb; j <= ub; j++)
-    {
-      if (iinf[j] > i)
-      {
-        iinf[j] = i;
-      }
-      if (isup[j] < i)
-      {
-	isup[j] = i; 
-
-      }
-    }
-  }
-
-	// determining extended mask edges for matrices 
+	const int rowmin = carte.rowmin();
+	const int rowmax = carte.rowmax();
+	const int colmin = carte.colmin();
+	const int colmax = carte.colmax();
+
+	// cout << "rowmin " << rowmin << endl;
+	// cout << "rowmax " << rowmax << endl;
+	// cout << "colmin " << colmin << endl;
+	// cout << "colmax " << colmax << endl;
+
+	IVECTOR iinftmp(colmin, colmax);
+	IVECTOR isuptmp(colmin, colmax);
+	iinftmp.initialize();
+	isuptmp.initialize();
+
+	// cout << "iinftmp " << iinftmp << endl;
+	// cout << "isuptmp " << isuptmp << endl;
+
+	for (int j = colmin; j <= colmax; j++)
+	{
+		int i = rowmin;
+		while (i <= rowmax && carte[i][j] == 0)
+		{
+			i++;
+		}
+		iinftmp[j] = i;
+	}
+
+	// cout << "iinftmp " << iinftmp << endl;
+
+	for (int j = colmax; j >= colmin; j--)
+	{
+		int i = rowmax;
+		while (i >= rowmin && carte[i][j] == 0)
+		{
+			i--;
+		}
+		isuptmp[j] = i;
+	}
+
+	// cout << "isuptmp" << isuptmp << endl;
+
+	jmin = colmin;
+	while (iinftmp(jmin) > isuptmp(jmin))
+		jmin++;
+
+	jmax = colmax;
+	while (iinftmp(jmax) > isuptmp(jmax))
+		jmax--;
+
+	for (int j = jmin; j <= jmax; j++)
+	{
+		if (iinftmp[j]>isuptmp[j]){
+			cerr << "Error: There is a vertical hole in the mask" << endl;
+			exit(2);
+		}
+	}
+
+	// cout << "jmin,jmax " << jmin << ' ' << jmax << endl;
+
+	for (int j = jmin; j <= jmax; j++)
+	{
+		int i = iinftmp[j];
+		if (!(carte[i - 1][j] == 0 && carte[i][j] != 0))
+		{
+			// cout << carte[i][j] << ' ' << carte[i + 1][j] << endl;
+			// cout << i << ' ' << j << endl;
+			// cout << __FILE__ << ':' << __LINE__ << endl;
+			exit(1);
+		}
+
+		i = isuptmp[j];
+		if (!(carte[i][j] != 0 && carte[i + 1][j] == 0))
+		{
+			// cout << carte[i][j] << ' ' << carte[i + 1][j] << endl;
+			// cout << i << ' ' << j << endl;
+			// cout << __FILE__ << ':' << __LINE__ << endl;
+			exit(1);
+		}
+	}
+
+	iinf.allocate(jmin, jmax);
+	isup.allocate(jmin, jmax);
+	iinf.initialize();
+	isup.initialize();
+
+	for (int j = jmin; j <= jmax; j++)
+	{
+		iinf[j] = iinftmp[j];
+		isup[j] = isuptmp[j];
+	}
+	// cout << "iinf " << iinf << endl;
+	// cout << "isup " << isup << endl;
+	// exit(1);
+	//=====================================
+
+	IVECTOR jinftmp(rowmin, rowmax);
+	IVECTOR jsuptmp(rowmin, rowmax);
+	jinftmp.initialize();
+	jsuptmp.initialize();
+	// cout << "jinftmp " << jinftmp << endl;
+	// cout << "jsuptmp " << jsuptmp << endl;
+
+	for (int i = rowmin; i <= rowmax; i++)
+	{
+		int j = colmin;
+		while (j <= colmax && carte[i][j] == 0)
+		{
+			j++;
+		}
+		jinftmp[i] = j;
+	}
+	// cout << "jinftmp " << jinftmp << endl;
+	// exit(1);
+
+	for (int i = rowmax; i >= rowmin; i--)
+	{
+		int j = colmax;
+		while (j >= colmin && carte[i][j] == 0)
+		{
+			j--;
+		}
+		jsuptmp[i] = j;
+	}
+	// cout << "jsuptmp " << jsuptmp << endl;
+	// exit(1);
+
+	imin = rowmin;
+	while (jinftmp(imin) > jsuptmp(imin))
+		imin++;
+	imax = rowmax;
+	while (jinftmp(imax) > jsuptmp(imax))
+		imax--;
+	// cout << imin << endl;
+	// cout << imax << endl;
+	// exit(1);
+
+	for (int i = imin; i <= imax; i++)
+	{
+		if (jinftmp[i]>jsuptmp[i]){
+			cerr << "Error: There is a horizontal hole in the mask" << endl;
+			exit(2);
+		}
+	}
+
+	jinf.allocate(imin, imax);
+	jsup.allocate(imin, imax);
+	for (int i = imin; i <= imax; i++)
+	{
+		jinf[i] = jinftmp[i];
+		jsup[i] = jsuptmp[i];
+	}
+	// cout << "jinf " << jinf << endl;
+	// cout << "jsup " << jsup << endl;
+
+	for (int j = jmin; j <= jmax; j++)
+	{
+		int lb = iinf[j];
+		int ub = isup[j];
+		for (int i = lb; i <= ub; i++)
+		{
+			if (jinf[i] > j)
+			{
+				jinf[i] = j;
+			}
+			if (jsup[i] < j)
+			{
+				jsup[i] = j;
+			}
+		}
+	}
+
+	// cout << "jinf " << jinf << endl;
+	// cout << "jsup" << jsup << endl;
+	// exit(1);
+
+	for (int i = imin; i <= imax; i++)
+	{
+		int lb = jinf[i];
+		int ub = jsup[i];
+		for (int j = lb; j <= ub; j++)
+		{
+			if (iinf[j] > i)
+			{
+				iinf[j] = i;
+			}
+			if (isup[j] < i)
+			{
+				isup[j] = i;
+			}
+		}
+	}
+
+	// determining extended mask edges for matrices
 	// used in computing finite differences on bounds
-  	imin1 = imin-1;
-	int jmin1 = jmin-1;
-	if (imin1 < rowmin) imin1 = rowmin;
-	if (jmin1 < colmin) jmin1 = colmin;
-	imax1 = imax+1;
-	int jmax1 = jmax+1;
-	if (imax1 > rowmax) imax1 = rowmax;
-	if (jmax1 > colmax) jmax1 = colmax;
-
-	jinf1.allocate(imin1,imax1); jinf1.initialize();
-	jsup1.allocate(imin1,imax1); jsup1.initialize();
-	
-	for (int j = colmin; j <= colmax; j++){
+	imin1 = imin - 1;
+	int jmin1 = jmin - 1;
+	if (imin1 < rowmin)
+		imin1 = rowmin;
+	if (jmin1 < colmin)
+		jmin1 = colmin;
+	imax1 = imax + 1;
+	int jmax1 = jmax + 1;
+	if (imax1 > rowmax)
+		imax1 = rowmax;
+	if (jmax1 > colmax)
+		jmax1 = colmax;
+
+	jinf1.allocate(imin1, imax1);
+	jinf1.initialize();
+	jsup1.allocate(imin1, imax1);
+	jsup1.initialize();
+
+	for (int j = colmin; j <= colmax; j++)
+	{
 		int i = rowmin;
-		while (i < rowmax && carte[i+1][j] == 0) i++;
+		while (i < rowmax && carte[i + 1][j] == 0)
+			i++;
 		iinftmp[j] = i;
 	}
 
-	for (int j = colmax; j >= colmin; j--){
+	for (int j = colmax; j >= colmin; j--)
+	{
 		int i = rowmax;
-		while (i > rowmin && carte[i-1][j] == 0) i--;
+		while (i > rowmin && carte[i - 1][j] == 0)
+			i--;
 		isuptmp[j] = i;
 	}
 
-	for (int i = imin1; i <= imax1; i++){
+	for (int i = imin1; i <= imax1; i++)
+	{
 		jinf1[i] = jinftmp[i];
 		jsup1[i] = jsuptmp[i];
 	}
 
-	for (int j = jmin1; j <= jmax1; j++) {
+	for (int j = jmin1; j <= jmax1; j++)
+	{
 		int lb = iinftmp[j];
 		int ub = isuptmp[j];
-		for (int i = lb; i <= ub; i++){
-			if (jinf1[i] > j-1){
-        			jinf1[i] = j-1;
-      			}
-      			if (jsup1[i] < j+1){
-				jsup1[i] = j+1; 
-      			}
-    		}
-  	}
+		for (int i = lb; i <= ub; i++)
+		{
+			if (jinf1[i] > j - 1)
+			{
+				jinf1[i] = j - 1;
+			}
+			if (jsup1[i] < j + 1)
+			{
+				jsup1[i] = j + 1;
+			}
+		}
+	}
 	/////////////////////////////////////////////////
 
-//cout << rowmin << " " << rowmax << " " << colmin << " " << colmax << endl;
-//cout << imin << " " << imax << " " << imin1 << " " << imax1 << endl;
-//cout << "jinf " << jinf << endl;
-//cout << "jsup" << jsup << endl;
-//cout << "jinf1 " << jinf1 << endl;
-//cout << "jsup1" << jsup1 << endl;
-//exit(1);
-
-/*
-  DMATRIX dmi(imin, imax, jinf, jsup);
-  dmi.initialize();
-  int count = 0;
-  for (int i = imin; i <= imax; i++)
-  {
-    const int nonzeros = jsup[i] - jinf[i] + 1;
-    count += jmax - nonzeros;
-cout << i << ' ' << nonzeros << ' ' << count << endl;
-    for (int j = jinf[i]; j <= jsup[i]; j++)
-    {
-      if (iinf[j] <= i && i <= isup[j])
-      {
-        dmi[i][j] += 1;
-      }
-    }
-  }
-//cout << dmi << endl;
-cout << imin << endl;
-cout << imax << endl;
-cout << jmin << endl;
-cout << jmax << endl;
-cout << count << endl;
-exit(1);
-  DMATRIX dmj(jmin, jmax, iinf, isup);
-  dmj.initialize();
-  for (int j = jmin; j <= jmax; j++)
-  {
-    for (int i = iinf[j]; i <= isup[j]; i++)
-    {
-      if (jinf[i] <= j && j <= jsup[i])
-      {
-        dmj[j][i] += 1;
-      }
-    }
-  }
-cout << dmj << endl;
-exit(1);
-*/
+	// cout << rowmin << " " << rowmax << " " << colmin << " " << colmax << endl;
+	// cout << imin << " " << imax << " " << imin1 << " " << imax1 << endl;
+	// cout << "jinf " << jinf << endl;
+	// cout << "jsup" << jsup << endl;
+	// cout << "jinf1 " << jinf1 << endl;
+	// cout << "jsup1" << jsup1 << endl;
+	// exit(1);
+
+	/*
+	  DMATRIX dmi(imin, imax, jinf, jsup);
+	  dmi.initialize();
+	  int count = 0;
+	  for (int i = imin; i <= imax; i++)
+	  {
+		const int nonzeros = jsup[i] - jinf[i] + 1;
+		count += jmax - nonzeros;
+	cout << i << ' ' << nonzeros << ' ' << count << endl;
+		for (int j = jinf[i]; j <= jsup[i]; j++)
+		{
+		  if (iinf[j] <= i && i <= isup[j])
+		  {
+			dmi[i][j] += 1;
+		  }
+		}
+	  }
+	//cout << dmi << endl;
+	cout << imin << endl;
+	cout << imax << endl;
+	cout << jmin << endl;
+	cout << jmax << endl;
+	cout << count << endl;
+	exit(1);
+	  DMATRIX dmj(jmin, jmax, iinf, isup);
+	  dmj.initialize();
+	  for (int j = jmin; j <= jmax; j++)
+	  {
+		for (int i = iinf[j]; i <= isup[j]; i++)
+		{
+		  if (jinf[i] <= j && j <= jsup[i])
+		  {
+			dmj[j][i] += 1;
+		  }
+		}
+	  }
+	cout << dmj << endl;
+	exit(1);
+	*/
 }
 
-void PMap::domain_type(const int nlon) {
+void PMap::domain_type(const int nlon)
+{
 	global = 0;
-	if (nlon==360) 
+	if (nlon == 360)
 		global = 1;
 
 	cout << "Is global domain? ";
-	if (global) cout << "YES" << endl;	
-	else cout << "NO" << endl;	
+	if (global)
+		cout << "YES" << endl;
+	else
+		cout << "NO" << endl;
 }
-
diff --git a/src/Map.h b/src/Map.h
index 4de7218..3d9bab1 100755
--- a/src/Map.h
+++ b/src/Map.h
@@ -25,6 +25,7 @@ class PMap
 
 public:
 	void lit_map(CParam &param);
+	void lit_tagmap(CParam &param, int tagpop);
 	void delete_map(const CParam &param) {/*DoesNothing*/}
 	void reg_indices(CParam &param);
 
@@ -35,8 +36,8 @@ class PMap
 
 	IMATRIX bord_cell;	// bord des cellules (voir Classe CBord plus loin)
 	IMATRIX nbl_bord_cell;	// border for the cells with all layers being exist
-	IMATRIX carte;		// carte (mettre des zéros quand terre ou ile)
-	DMATRIX itopo;		// carte (mettre des zéros quand terre ou ile)
+	IMATRIX carte;		// carte (mettre des z�ros quand terre ou ile)
+	DMATRIX itopo;		// carte (mettre des z�ros quand terre ou ile)
 	IMATRIX maskEEZ;	// carte des EEZ avec un code par EEZ
 	IMATRIX maskMPA;	// carte des MPA avec un code par MPA
 
@@ -62,9 +63,9 @@ class PMap
 	ivector regjmax;
 
 private:
-	void definit_cell_bords(const int nti, const int ntj);	// définition des bords des cellules
-	void definit_lim_infsup(const int nti, const int ntj);	// définition des limites inférieures 
-								// et supérieures de la matrice zone pour optimiser le calcul
+	void definit_cell_bords(const int nti, const int ntj);	// d�finition des bords des cellules
+	void definit_lim_infsup(const int nti, const int ntj);	// d�finition des limites inf�rieures 
+								// et sup�rieures de la matrice zone pour optimiser le calcul
 	void domain_type(const int nlon); 
 };
 
diff --git a/src/Param.h b/src/Param.h
index 68f08ad..cb75f0a 100755
--- a/src/Param.h
+++ b/src/Param.h
@@ -140,6 +140,7 @@ class CParam
 	string  str_dir_init;
 	string  str_dir_fisheries;
 	string  str_dir_tags;
+	string 	str_dir_tagmasks;
 	string  strfile_pp;
 	string  strfile_sst;
 	string  strfile_vld;
@@ -157,6 +158,7 @@ class CParam
 	vector<string> strfile_umc,strfile_vmc,strfile_tmc, strfile_oxymc;		// array of string [nb_layer]
 	string strdir_output;
 	string strout_tags;
+	int use_tag_masks; // whether to use individual masks for tag cohorts
 	int write_all_cohorts_dym;
 	int write_all_fisheries_dym;
 	
@@ -279,6 +281,7 @@ class CParam
 	vector<string>  file_catch_data;
 	vector<string>  file_frq_data;
 	vector<string>  file_tag_data;
+	vector<string>  file_tag_masks;
 	int nb_catch_files, nb_frq_files, nb_tag_files;
 	int tag_gauss_kernel_on;
 	float dx_tags, dy_tags; 	// setup of the grid to aggregate tagging data
diff --git a/src/SeapodymCoupled.h b/src/SeapodymCoupled.h
index a006652..632a250 100755
--- a/src/SeapodymCoupled.h
+++ b/src/SeapodymCoupled.h
@@ -162,8 +162,9 @@ friend class tag_release;
 	//void Age_class_survivals(dvar_matrix& N_a, const dmatrix N_a_1, const int nt_a, const int nt_a_1);
 	///void Survival(dvar_matrix& N_a, dvar_matrix& N_a_1, dvar_matrix& Mort_a, dvar_matrix& Mort_a_1, const int a, const int sp, const bool adult);
 	void Survival(dvar_matrix& N_a, dvar_matrix& N_a_1, const int a, const int sp);
-	void Ageing(dvar_matrix& N_a, dvar_matrix& N_a_1);
-        void AgePlus(dvar_matrix& N_a, dvar_matrix& N_a_1);
+	void Survival_tagpop(dvar_matrix& N_a, dvar_matrix& N_a_1, const int a, const int sp, const int tagpop);
+	void Ageing(dvar_matrix& N_a, dvar_matrix& N_a_1, PMap& map);
+    void AgePlus(dvar_matrix& N_a, dvar_matrix& N_a_1, PMap& map);
 
 	void Spawning(dvar_matrix& J, dvar_matrix& Hs, dvar_matrix& N_a, const int jday, const int sp, const int pop_built, const int t_count);
 	void spawning_spinup_comp(dmatrix& J, dmatrix& Hs, dmatrix& SST, double nb_recruitment, double Tmin);
diff --git a/src/SeapodymCoupled_EditRunCoupled.cpp b/src/SeapodymCoupled_EditRunCoupled.cpp
index 719de17..b443c33 100755
--- a/src/SeapodymCoupled_EditRunCoupled.cpp
+++ b/src/SeapodymCoupled_EditRunCoupled.cpp
@@ -39,6 +39,15 @@ int SeapodymCoupled::EditRunCoupled(const char* parfile)
 				.5*param->sp_unit_cohort[sp][a-1]+.5*param->sp_unit_cohort[sp][a];
 	}
 	map.lit_map(*param);
+	nb_tagpops = param->nb_tag_files;
+
+	if (param->tag_like[0] && param->use_tag_masks){
+		for (int tagpop=0; tagpop<nb_tagpops; tagpop++){
+			PMap tagmap;
+			tagmap.lit_tagmap(*param, tagpop);
+			tagmaps.push_back(tagmap);
+		}
+	}
 	pop.InitCalPop(*param, map);
 	//Notice, here param variables will be rewritten, it's better to add the test to check their identity
 	//rw.rbin_headpar(param->strfile_pp, param->nlong, param->nlat, param->nlevel);
diff --git a/src/SeapodymCoupled_Funcs.cpp b/src/SeapodymCoupled_Funcs.cpp
index fa62159..d63ac07 100644
--- a/src/SeapodymCoupled_Funcs.cpp
+++ b/src/SeapodymCoupled_Funcs.cpp
@@ -254,12 +254,23 @@ void SeapodymCoupled::CalcSums()
 					if (param->tag_like[sp]){
 						for (int aa=a0_adult[sp]; aa<aN_adult[sp]; aa++){
 							double total_tags = 0.0;
-							for (int pop=1; pop<=param->nb_tag_files;pop++)
-								total_tags += value(mat.dvarDensity(pop,aa,i,j));
+							for (int pop=1; pop<=param->nb_tag_files;pop++){
+
+								int add_to_total_tags = 0;
+								if (!param->use_tag_masks){
+									add_to_total_tags = 1;
+								}else if ((i > tagmaps[pop-1].imin1 && i < tagmaps[pop-1].imax1 && j > tagmaps[pop-1].jinf[i] && j < tagmaps[pop-1].jsup[i]))
+								{
+									add_to_total_tags = 1;
+								}
+								if (add_to_total_tags){
+									total_tags += value(mat.dvarDensity(pop,aa,i,j));
+
+									if (pop==1 && value(mat.dvarDensity(pop,aa,i,j))<0)  cerr << "NEGATIVE BIOMASS for " << aa << " "  << value(mat.dvarDensity(pop,aa,i,j)) << endl;
+								}
+							}
 
 							mat.recruit[0][i][j] += total_tags; 
-							if (value(mat.dvarDensity(1,aa,i,j))<0) cerr << "NEGATIVE BIOMASS for " << aa << " " 
-												     << value(mat.dvarDensity(1,aa,i,j)) << endl;
 						}
 					}
 					//mat.sum_B_recruit[sp] += value(mat.dvarDensity[sp][age_recruits][i][j])*W_mt(age_recruits)*lat_corrected_area;
diff --git a/src/SeapodymCoupled_OnRunCoupled.cpp b/src/SeapodymCoupled_OnRunCoupled.cpp
index 530ce82..14fbc50 100755
--- a/src/SeapodymCoupled_OnRunCoupled.cpp
+++ b/src/SeapodymCoupled_OnRunCoupled.cpp
@@ -432,10 +432,18 @@ double SeapodymCoupled::OnRunCoupled(dvar_vector x, const bool writeoutputfiles)
 							if (!tagpop_age_solve(spop-1,t_count-1,age))
 								continue;
 						}
-						pop.Precalrec_Calrec_adult(map,mat,*param,rw,
-							mat.dvarDensity[spop][age],Mortality,
-							tcur,fishing,age,sp,year,month,jday,
-							step_fishery_count,mortality_off);//checked 20150210
+						
+						PMap* map_ptr=nullptr;
+						if (spop>0 && param->use_tag_masks){
+							if ((spop-1) < static_cast<int>(tagmaps.size()))
+								map_ptr = &tagmaps[spop-1]; 
+						}else{
+							map_ptr = &map;
+						}
+						pop.Precalrec_Calrec_adult(*map_ptr,mat,*param,rw,
+						mat.dvarDensity[spop][age],Mortality,
+						tcur,fishing,age,sp,year,month,jday,
+						step_fishery_count,mortality_off);//checked 20150210
 					}
 					if (sum(param->mask_fishery_sp) && !tags_only) fishing = true;
 				}
@@ -478,7 +486,11 @@ double SeapodymCoupled::OnRunCoupled(dvar_vector x, const bool writeoutputfiles)
 							}
 						}
 						//6.2. Ageing of population density
-						Survival(mat.dvarDensity[spop][a], mat.dvarDensity[spop][a-1] , a, sp);
+						if (spop>0 && param->use_tag_masks){
+							Survival_tagpop(mat.dvarDensity[spop][a], mat.dvarDensity[spop][a-1] , a, sp, spop-1);
+						}else{
+							Survival(mat.dvarDensity[spop][a], mat.dvarDensity[spop][a-1] , a, sp);
+						}
 					}
 				}
 				nt_dtau=0;
@@ -551,7 +563,6 @@ double SeapodymCoupled::OnRunCoupled(dvar_vector x, const bool writeoutputfiles)
 
 	} // end of simulation loop
 
-
 	if (writeoutputfiles) {SaveDistributions(year, month);
 		cout << "total catch in simulation (optimization): " << SUM_CATCH << endl;
 	}
diff --git a/src/SeapodymCoupled_OnRunFirstStep.cpp b/src/SeapodymCoupled_OnRunFirstStep.cpp
index d9cacba..ac5a19b 100755
--- a/src/SeapodymCoupled_OnRunFirstStep.cpp
+++ b/src/SeapodymCoupled_OnRunFirstStep.cpp
@@ -24,7 +24,11 @@ void SeapodymCoupled::OnRunFirstStep()
 	mat.createMatForage(map, nb_forage, t0, nbt, nbi, nbj);
 
 	int nb_pops = nb_species*(1+param->nb_tag_files);
-	mat.CreateMatSpecies(map,t0, nbt, nbi, nbj, nb_pops, a0_adult, param->sp_nb_cohorts);
+	if (param->use_tag_masks){
+		mat.CreateMatSpecies(map, tagmaps, t0, nbt, nbi, nbj, nb_pops, a0_adult, param->sp_nb_cohorts);
+	}else{
+		mat.CreateMatSpecies(map,t0, nbt, nbi, nbj, nb_pops, a0_adult, param->sp_nb_cohorts);
+	}
 	mat.CreateMatHabitat(map,nb_species,nb_forage,nb_layer,max(param->sp_nb_cohorts),t0, nbt,nbi,nbj,a0_adult,param->sp_nb_cohorts,param->age_compute_habitat);
 	past_month=0;
 	past_qtr=0;
@@ -78,7 +82,6 @@ void SeapodymCoupled::OnRunFirstStep()
 	}
 	func.allocate_dvmatr(map.imin,map.imax,map.jinf,map.jsup);
 	//TAG data reading and allocation section
-	nb_tagpops = param->nb_tag_files;
 
 	if (param->tag_like[0]){
 		nb_rel.allocate(0,nb_tagpops-1);
diff --git a/src/SeapodymCoupled_ReadTags.cpp b/src/SeapodymCoupled_ReadTags.cpp
index 4ca3960..f99880e 100755
--- a/src/SeapodymCoupled_ReadTags.cpp
+++ b/src/SeapodymCoupled_ReadTags.cpp
@@ -141,6 +141,13 @@ void SeapodymCoupled::ReadTaggingData(imatrix& nb_rel, ivector& t_count_rec)
 
 	int nbtot_tags_files = 0;
 	for (int p=0; p<nb_tagpops; p++){
+		PMap* map_ptr=nullptr;
+		if (param->use_tag_masks){
+			if (p < static_cast<int>(tagmaps.size()))
+				map_ptr = &tagmaps[p];
+		}else{
+			map_ptr = &map;
+		}
 
 		string file_in = param->file_tag_data[p];
 		//date of all recaptures in the cohort will be read from the name of the file
@@ -208,7 +215,7 @@ void SeapodymCoupled::ReadTaggingData(imatrix& nb_rel, ivector& t_count_rec)
 				j_mod = param->lattoj(lat_rec);	
 				if (lon_rec>xlon[0] && lon_rec<=xlon[nx_obs] 
 				    && lat_rec<=ylat[0] && lat_rec>ylat[ny_obs]){
-					if (map.carte(i_mod,j_mod))
+					if (map_ptr->carte(i_mod,j_mod))
 						rec_added = 1;
 				}
 
@@ -217,7 +224,7 @@ void SeapodymCoupled::ReadTaggingData(imatrix& nb_rel, ivector& t_count_rec)
 				i_mod = param->lontoi(lon);
 				j_mod = param->lattoj(lat);
 				if (i_mod>=map.imin && i_mod<=map.imax && j_mod>=map.jmin && j_mod<=map.jmax){
-					if (map.carte(i_mod,j_mod))
+					if (map_ptr->carte(i_mod,j_mod))
 						rel_added = 1;
 					//attn: generic way is to search over 8 surrounding cells
 					//if (!map.carte(i_mod,j_mod) && map.carte(i_mod,j_mod-1)) { 
@@ -244,7 +251,7 @@ void SeapodymCoupled::ReadTaggingData(imatrix& nb_rel, ivector& t_count_rec)
 							cout << "WARNING: tag recapture is on land: " << 
 								p << " " << id << " " << yy_rec << " "<< mm_rec 
 								<< " " << lat_rec << " " << lon_rec << "; mask: " 
-								<< map.carte(i,j)<< endl;
+								<< map_ptr->carte(i,j)<< endl;
 						} else 
 							cout<< "WARNING: tag recapture out of observational space: " << 
 								p << " " << id << " " << yy_rec << " "<< mm_rec 
@@ -258,7 +265,7 @@ void SeapodymCoupled::ReadTaggingData(imatrix& nb_rel, ivector& t_count_rec)
 						cout<< "WARNING: tag release out of model domain: " << 
 							p << " " << id << " " << yy << " " << mm << " "
 							<< lat << " " << lon << "; mask: " 
-							<< map.carte(i_mod,j_mod) << endl;
+							<< map_ptr->carte(i_mod,j_mod) << endl;
 					continue;
 				}
 				//One more check might be necessary: discrepancy between length(age) at recapture
@@ -325,7 +332,7 @@ void SeapodymCoupled::ReadTaggingData(imatrix& nb_rel, ivector& t_count_rec)
 							int jto   = jfrom + yr_tags;
 							for (int i=ifrom; i<ito; i++){
 								for (int j=jfrom; j<jto; j++){
-									if (map.carte[i][j]){
+									if (map_ptr->carte[i][j]){
 										rec_obs(p,ii,jj) += gkernel(i,j);
 										tlib_obs(p,ii,jj) += gkernel(i,j)*dtlib;
 									}
diff --git a/src/SeapodymDocConsole.h b/src/SeapodymDocConsole.h
index 58bc39e..68c11f8 100755
--- a/src/SeapodymDocConsole.h
+++ b/src/SeapodymDocConsole.h
@@ -26,6 +26,7 @@ class SeapodymDocConsole
 	VarParamCoupled* param;	// objet param derive de la classe CParam 
 	VarMatrices mat;	// objet mat derive de la classe CMatrices
 	PMap map;		// objet map derive de la classe CMap
+	vector<PMap> tagmaps;		// tableau de maps for tag populations
 	VarSimtunaFunc func;	// objet func derive de la classe CSimtunafunc
 	CNumfunc nfunc; 	// objet nfunc derive de la classe CNumfunc
 	CSaveTimeArea save;	// objet save derive de la classe CSaveTimeArea
diff --git a/src/VarMatrices.h b/src/VarMatrices.h
index d69b29e..eb77c7a 100755
--- a/src/VarMatrices.h
+++ b/src/VarMatrices.h
@@ -75,9 +75,9 @@ class VarMatrices : public CMatrices
 
 			for (int age = 0; age < agemax_sp; age++) {
 				dvarDensity(sp, age).allocate(map.imin1, map.imax1, map.jinf1, map.jsup1);
-				dvarDensity(sp, age).initialize();
 			}
 		}
+
 /*
 		const int age_max_1 = 48-1;
 		dvarDensity_age.allocate(0,age_max-1);
@@ -89,6 +89,27 @@ class VarMatrices : public CMatrices
 
 	}
 
+	void CreateMatSpecies(PMap& map, vector<PMap> tagmaps, int t0, int nbt, int nbi, int nbj, int nb_species, const ivector a0_adult, const ivector& sp_nb_cohorts) {
+		//CMatrices::createMatSpecies(map, t0, nbt, nbi, nbj, nb_species, a0_adult, sp_nb_cohorts);
+		CMatrices::createMatSpecies(map, t0, nbt, nbi, nbj, 1, a0_adult, sp_nb_cohorts);
+		dvarDensity.allocate(0, nb_species - 1);
+		for (int sp = 0; sp < nb_species; sp++) {
+			//const int agemax_sp = sp_nb_cohorts(sp);
+			const int agemax_sp = sp_nb_cohorts(0);
+			dvarDensity(sp).allocate(0, agemax_sp - 1);
+			dvarDensity(sp).initialize();
+
+			for (int age = 0; age < agemax_sp; age++) {
+				if (sp==0){
+					dvarDensity(sp, age).allocate(map.imin1, map.imax1, map.jinf1, map.jsup1);
+				}else{
+					dvarDensity(sp, age).allocate(tagmaps[sp-1].imin1, tagmaps[sp-1].imax1, tagmaps[sp-1].jinf1, tagmaps[sp-1].jsup1);
+				}
+				dvarDensity(sp, age).initialize();
+			}
+		}
+	}
+
 	void CreateMatCatch(PMap& map, int nbi, int nbj, int nb_species, const IVECTOR& nb_fleet, const ivector a0_adult, const IVECTOR& nb_cohorts, const IVECTOR& nb_region) {
 		CMatrices::createMatCatch(map, nbi, nbj, nb_species, nb_fleet, a0_adult, nb_cohorts, nb_region);
 
diff --git a/src/VarParamCoupled.cpp b/src/VarParamCoupled.cpp
index d722388..0287dbd 100755
--- a/src/VarParamCoupled.cpp
+++ b/src/VarParamCoupled.cpp
@@ -115,7 +115,8 @@ bool VarParamCoupled::read(const string& parfile)
 	str_dir_tags = str_dir;
 	if (!doc.get("/strdir_tags", "value").empty())
 		str_dir_tags = doc.get("/strdir_tags", "value");
-
+	if (!doc.get("/strdir_tagmasks", "value").empty())
+		str_dir_tagmasks = doc.get("/strdir_tagmasks", "value");
 	strfile_pp = str_dir + doc.get("/strfile_pp", "value");
 
 	use_sst = 0;
@@ -990,6 +991,11 @@ bool VarParamCoupled::read(const string& parfile)
 				lonmin_tags = 180; lonmax_tags = 290;latmin_tags = -50;latmax_tags = 3;
 			}
 			nb_tag_files = 0;
+			if (!doc.get("/use_tag_masks", "flag").empty()){
+				use_tag_masks=doc.getInteger("/use_tag_masks", "flag");
+			}else{
+				use_tag_masks=0;
+			}
 			if (tag_like(sp)){
 				if (!doc.get("/tags_grid").empty()){
 					dx_tags = doc.getDouble(string("/tags_grid/reso"),"dx");
@@ -1015,7 +1021,15 @@ bool VarParamCoupled::read(const string& parfile)
 					for (int nf=0; nf<nb_tag_files; nf++){
 						std::ostringstream ostr;
 						ostr << "file" << nf+1;
-						file_tag_data.push_back(str_dir_tags+doc.get("/file_tag_data/"+sp_name[sp],ostr.str()));
+						string file = doc.get("/file_tag_data/"+sp_name[sp],ostr.str());
+						file_tag_data.push_back(str_dir_tags+file);
+						if (use_tag_masks){
+							if (!str_dir_tagmasks.empty()){
+								file_tag_masks.push_back(str_dir_tagmasks+file);
+							}else{
+								cout << "WARNING: no <strdir_tagmasks> field in parameter file" << endl; 
+							}
+						}
 					}
 				} else {
 					cout << "WARNING: TAG DATA ARE ABSENT" << endl; 
diff --git a/src/calrec_adre.cpp b/src/calrec_adre.cpp
index 83a83fc..7207649 100755
--- a/src/calrec_adre.cpp
+++ b/src/calrec_adre.cpp
@@ -35,7 +35,8 @@ void CCalpop::calrec1(const PMap& map, dvar_matrix& uu, const dmatrix& mortality
 
 					rhs[i]=-d[i][j]*value(uu(i,j-1)) + (2*iterationNumber-e[i][j])*value(uu(i,j)) - f[i][j]*value(uu(i,j+1));
 
-                         	} else {
+                } else {
+					cerr << "Error: cf file " << __FILE__ << ", line " << __LINE__ << ": itr = " << itr << ", i=" << i << ",j=" << j <<  endl;
 					cout << __LINE__ << endl; exit(1);
 					//rhs[i]=(2*iterationNumber)*value(uu(i,j));
 				}
diff --git a/src/dv_survival.cpp b/src/dv_survival.cpp
index b9e672f..5c1c69e 100755
--- a/src/dv_survival.cpp
+++ b/src/dv_survival.cpp
@@ -16,19 +16,27 @@ unsigned long int restore_long_int_value(void);
 void SeapodymCoupled::Survival(dvar_matrix& N_a, dvar_matrix& N_a_1, const int a, const int sp)
 {
 	int nb_cohorts = param->sp_nb_cohorts[sp];
+	PMap* map_ptr=&map;
 
 	if (a<nb_cohorts-1)
-
-		Ageing(N_a, N_a_1);
+		Ageing(N_a, N_a_1, *map_ptr);
 	else 
+		AgePlus(N_a, N_a_1, *map_ptr);
+} 
 
-		AgePlus(N_a, N_a_1);
-	
-	
+void SeapodymCoupled::Survival_tagpop(dvar_matrix& N_a, dvar_matrix& N_a_1, const int a, const int sp, const int tagpop)
+{
+	int nb_cohorts = param->sp_nb_cohorts[sp];
+	PMap* map_ptr = &tagmaps[tagpop];
 
+	if (a<nb_cohorts-1)
+		Ageing(N_a, N_a_1, *map_ptr);
+	else 
+		AgePlus(N_a, N_a_1, *map_ptr);
 } 
 
-void SeapodymCoupled::Ageing(dvar_matrix& N_a, dvar_matrix& N_a_1)
+
+void SeapodymCoupled::Ageing(dvar_matrix& N_a, dvar_matrix& N_a_1, PMap& map)
 {
 	for (int i = map.imin; i <= map.imax; i++){	
 		const int jmin = map.jinf[i];
@@ -51,7 +59,7 @@ void SeapodymCoupled::Ageing(dvar_matrix& N_a, dvar_matrix& N_a_1)
 
 }
 
-void SeapodymCoupled::AgePlus(dvar_matrix& N_a, dvar_matrix& N_a_1)
+void SeapodymCoupled::AgePlus(dvar_matrix& N_a, dvar_matrix& N_a_1, PMap& map)
 {
 	N_a.save_dvar_matrix_position();
 	save_identifier_string2((char*)"before_ageplus");
diff --git a/src/fd_survival.cpp b/src/fd_survival.cpp
index 0ca86ad..c2a2c40 100755
--- a/src/fd_survival.cpp
+++ b/src/fd_survival.cpp
@@ -4,17 +4,26 @@
 ///Forward main and forward functions called in simulation mode only for: 
 ///ageing term discretisation either for all 0..n_a-1 or for A+ age class.
 
-void SeapodymCoupled::Survival(dvar_matrix& N_a, dvar_matrix& N_a_1, const int a, const int sp)
+oid SeapodymCoupled::Survival(dvar_matrix& N_a, dvar_matrix& N_a_1, const int a, const int sp)
 {
 	int nb_cohorts = param->sp_nb_cohorts[sp];
+	PMap* map_ptr=&map;
 
 	if (a<nb_cohorts-1)
-
-		Ageing(N_a, N_a_1);
+		Ageing(N_a, N_a_1, *map_ptr);
 	else 
+		AgePlus(N_a, N_a_1, *map_ptr);
+} 
 
-		AgePlus(N_a, N_a_1);
+void SeapodymCoupled::Survival_tagpop(dvar_matrix& N_a, dvar_matrix& N_a_1, const int a, const int sp, const int tagpop)
+{
+	int nb_cohorts = param->sp_nb_cohorts[sp];
+	PMap* map_ptr = &tagmaps[tagpop];
 
+	if (a<nb_cohorts-1)
+		Ageing(N_a, N_a_1, *map_ptr);
+	else 
+		AgePlus(N_a, N_a_1, *map_ptr);
 } 
 
 void SeapodymCoupled::Ageing(dvar_matrix& N_a, dvar_matrix& N_a_1)
diff --git a/src/like.cpp b/src/like.cpp
index 78cce0f..f2960a2 100755
--- a/src/like.cpp
+++ b/src/like.cpp
@@ -143,13 +143,20 @@ double SeapodymCoupled::get_tag_like(dvariable& likelihood, bool writeoutputs)
 		//int nb_obs = 0;
 		if (t_count==t_count_rec(p)){
 
-			const int imin = map.imin; 
-			const int imax = map.imax; 
+			PMap* map_ptr=nullptr;
+			if (param->use_tag_masks){
+				if (p < static_cast<int>(tagmaps.size()))
+					map_ptr = &tagmaps[p]; 
+			}else{
+				map_ptr = &map;
+			}
+			const int imin = map_ptr->imin; 
+			const int imax = map_ptr->imax; 
 			for (int i = imin; i <= imax; i++){
-				const int jmin = map.jinf[i];
-				const int jmax = map.jsup[i];
+				const int jmin = map_ptr->jinf[i];
+				const int jmax = map_ptr->jsup[i];
 				for (int j = jmin ; j <= jmax; j++){
-					if (map.carte[i][j]){
+					if (map_ptr->carte[i][j]){
 										
 						double xx = param->itolon(i);
 						double yy = param->jtolat(j);
@@ -166,10 +173,10 @@ double SeapodymCoupled::get_tag_like(dvariable& likelihood, bool writeoutputs)
 				}
 				//temporally write recaptures to catch matrix (comment in CalcSums)
 				for (int i = imin; i <= imax; i++){
-					const int jmin = map.jinf[i];
-					const int jmax = map.jsup[i];
+					const int jmin = map_ptr->jinf[i];
+					const int jmax = map_ptr->jsup[i];
 					for (int j = jmin ; j <= jmax; j++){
-						if (map.carte[i][j]){
+						if (map_ptr->carte[i][j]){
 							double xx = param->itolon(i);
 							double yy = param->jtolat(j);
 							if (writeoutputs){
@@ -239,7 +246,7 @@ double SeapodymCoupled::get_tag_like(dvariable& likelihood, bool writeoutputs)
 					if (month>9) ostr << month <<15;
 					else
 						ostr << 0 << month <<15;
-					string file_out = "./" + param->sp_name[0] + "_tags_pred_"  + ostr.str() + ".txt";
+					string file_out = param->strout_tags + param->sp_name[0] + "_tags_pred_"  + ostr.str() + ".txt";
 					wtxt.open(file_out.c_str(), ios::out);
 					
 					if (wtxt){
@@ -249,7 +256,7 @@ double SeapodymCoupled::get_tag_like(dvariable& likelihood, bool writeoutputs)
 					}
 					wtxt.close();
 			
-					file_out = "./" + param->sp_name[0] + "_tags_obs_"  + ostr.str() + ".txt";
+					file_out = param->strout_tags + param->sp_name[0] + "_tags_obs_"  + ostr.str() + ".txt";
 					//file_out = "./" + param->sp_name[0] + "_releases_obs_"  + ostr.str() + ".txt";
 					wtxt.open(file_out.c_str(), ios::out);
 					if (wtxt){
diff --git a/template_parfile.xml b/template_parfile.xml
new file mode 100644
index 0000000..ac2eeec
--- /dev/null
+++ b/template_parfile.xml
@@ -0,0 +1,712 @@
+<?xml version="1.0"?>
+<!-- check http://www.w3.org/TR/2004/REC-xml-20040204 for current xml specifications -->
+<par version="1.0">
+
+<!--==================== -->
+<!--SIMULATION PARAMETERS-->
+<!--==================== -->
+<latitudeMin value="-54.5"/> 
+<latitudeMax value="65.5"/>
+<longitudeMin value="88.5"/>
+<longitudeMax value="290.5"/>
+<space_reso value="2"/>
+<nb_layer value="3"/>
+<deltaX value="120"/>
+<deltaY value="120"/>
+<deltaT value="30"/>
+<iterationNumber value="15"/>
+<tuna_spinup value="0"/>
+<save_first_date year="1978" month="1"/>
+<save_last_date year="2010" month="12"/>
+<nb_yr_forecast value="0"/>
+<nb_step_to_skip value="48"/>
+
+<!--==================== -->
+<!--LIKELIHOOD PARAMETERS--> 
+<!--==================== -->
+
+<total_likelihood value="1057384.045513833"/>
+
+<tag_likelihood_only value="1"/>
+
+<!--type of data in the likelihood: 0 - cpue, 1 - catch -->
+<like_c_cpue value="1"/>
+
+<!--likelihood function: 1 - concentrated, 2 - log-normal, 3 - poisson, 4 - negative binomial, 5 - zero-inflated negbin, 6 - truncated poisson -->
+<!--can be different for each fishery (see below the fishery codes) -->
+<likelihood_types>
+	<skj>3 1 3 3 3 1 1 1 1 1 1 1 1 1 3</skj>
+</likelihood_types>
+
+<nb_fishery value="15"/>
+<list_fishery_name>P1 P21 P22 P23 P3 S4 S5 S6 S7 S10 S11 S12 S13 S14 P15</list_fishery_name> 
+
+<!--1 - C in metric tones; 0 - C in numbers -->
+<fishery_catch_units>1 1 1 1 1 1 1 1 1 1 1 1 1 1 1</fishery_catch_units> 
+
+<degrade_fishery_reso_deg value="5"/>
+
+<!--Mask for fishery data usage: 0 - do not take into account-->
+<mask_fishery_sp>
+	<skj>1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1</skj>
+	<!--<skj>0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0</skj>-->
+</mask_fishery_sp>
+
+<mask_fishery_no_effort>
+        <skj>0 1 0 0 0 1 1 1 1 1 1 1 1 1 0</skj>
+</mask_fishery_no_effort>
+
+<mask_fishery_likelihood>
+	<!--<skj>1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1</skj>-->
+	<!--<skj>1 0 1 1 1 1 1 0 1 1 1 1 1 1 1</skj>-->
+	<skj>0 0 0 0 0 0 0 0 0 0 0 0 0 0 0</skj>
+</mask_fishery_likelihood>
+
+<likelihood_type value="1"/>
+
+<!--length frequencies likelihood: 0 - do not take into account, 1 - account--> 
+<frq_likelihood skj="0"/>
+
+<frq_likelihood_weight skj="0.07"/>
+
+<!--TAG likelihood: 0 - do not take into account, 1 - account-->
+<tag_likelihood skj="1"/>
+
+<!--STOCK likelihood: 0 - do not take into account, 1 - account-->
+<stock_likelihood skj="1"/>
+<!-- The values entered for mean_stock_obs will be used in stock likelihood
+     These can be the total abundance estimated by stock assessment models 
+     within the rectangular region specified -->
+<mean_stock_obs>
+     <!-- 4850 is MFCL_estimate for 2011 in WCPO -->
+    <skj value="8000" lgmin="120" lgmax="290" ltmin="-20" ltmax="45"/>
+</mean_stock_obs>
+
+<!--================= -->
+<!--MODEL INPUT-OUTPUT-->
+<!--================= -->
+<!-- DataFiles -->
+
+<!-- Working Directory -->
+<strdir value="/home/lucasb/SEAPODYM/data/forcings/run-interim_2x30d_po/"/>
+
+<str_file_mask value="po_interim_2deg_mask_short_noIO.txt"/>
+<!--<str_file_topo value="po_interim_2deg_topo_noeff.txt"/>-->
+<str_file_topo value="po_interim_2deg_topo_100m_modif.txt"/>
+
+<strdir_init value="init/1979/skj/ref_interim/Fref/"/>
+<strdir_forage value="forage/"/>
+
+<strdir_fisheries value="/home/lucasb/SEAPODYM/data/fisheries/skj/ds2013/"/>
+
+<file_catch_data>
+       <nb_files skj="1"/>  <skj file1="skj_catch_17f_PO_config2.txt"/>
+
+</file_catch_data>
+
+<flex_regstruc value="1"/>
+<file_frq_data>
+<nb_files skj="2"/>  <skj file1="skj_LF_11f_WCPO_L8mod.txt" file2="skj_LF_5f_PS_EPO_filtered.txt"/>
+</file_frq_data>
+
+
+		<!--<strdir_tags value="./ds2016/"/>-->
+<strdir_tags value="/home/lucasb/SEAPODYM/data/tagging/skj/ds2016/"/>
+
+<tags_grid>
+        <reso dx="6" dy="6"/>
+        <longitude east="110" west="280"/>
+        <latitude south="-16" north="44"/>
+</tags_grid>
+
+<file_tag_data>
+		<nb_files skj="152"/>
+		<!--<skj file1="skj-pop-tags_20090115.txt" file2="skj-pop-tags_20090215.txt" file3="skj-pop-tags_20090315.txt" file4="skj-pop-tags_20090415.txt" file5="skj-pop-tags_20090515.txt" file6="skj-pop-tags_20090615.txt" file7="skj-pop-tags_20090815.txt" file8="skj-pop-tags_20091015.txt" file9="skj-pop-tags_20091115.txt" file10="skj-pop-tags_20091215.txt" file11="skj-pop-tags_20100115.txt" file12="skj-pop-tags_20100215.txt" file13="skj-pop-tags_20100315.txt" file14="skj-pop-tags_20100415.txt" file15="skj-pop-tags_20100515.txt" file16="skj-pop-tags_20100615.txt" file17="skj-pop-tags_20100715.txt" file18="skj-pop-tags_20100815.txt" file19="skj-pop-tags_20100915.txt" file20="skj-pop-tags_20101015.txt" file21="skj-pop-tags_20101115.txt" file22="skj-pop-tags_20101215.txt"/>-->
+        <skj file1="skj-pop-tags_19780115.txt" file2="skj-pop-tags_19780215.txt" file3="skj-pop-tags_19780315.txt" file4="skj-pop-tags_19780415.txt" file5="skj-pop-tags_19780515.txt" file6="skj-pop-tags_19780615.txt" file7="skj-pop-tags_19780715.txt" file8="skj-pop-tags_19780815.txt" file9="skj-pop-tags_19780915.txt" file10="skj-pop-tags_19781015.txt" file11="skj-pop-tags_19781115.txt" file12="skj-pop-tags_19781215.txt" file13="skj-pop-tags_19790115.txt" file14="skj-pop-tags_19790215.txt" file15="skj-pop-tags_19790315.txt" file16="skj-pop-tags_19790415.txt" file17="skj-pop-tags_19790515.txt" file18="skj-pop-tags_19790615.txt" file19="skj-pop-tags_19790715.txt" file20="skj-pop-tags_19790815.txt" file21="skj-pop-tags_19790915.txt" file22="skj-pop-tags_19791015.txt" file23="skj-pop-tags_19791115.txt" file24="skj-pop-tags_19791215.txt" file25="skj-pop-tags_19800115.txt" file26="skj-pop-tags_19800215.txt" file27="skj-pop-tags_19800315.txt" file28="skj-pop-tags_19800415.txt" file29="skj-pop-tags_19800515.txt" file30="skj-pop-tags_19800615.txt" file31="skj-pop-tags_19800715.txt" file32="skj-pop-tags_19800815.txt" file33="skj-pop-tags_19800915.txt" file34="skj-pop-tags_19801015.txt" file35="skj-pop-tags_19801115.txt" file36="skj-pop-tags_19801215.txt" file37="skj-pop-tags_19810115.txt" file38="skj-pop-tags_19810215.txt" file39="skj-pop-tags_19810315.txt" file40="skj-pop-tags_19810415.txt" file41="skj-pop-tags_19810515.txt" file42="skj-pop-tags_19810615.txt" file43="skj-pop-tags_19810715.txt" file44="skj-pop-tags_19810815.txt" file45="skj-pop-tags_19810915.txt" file46="skj-pop-tags_19811215.txt" file47="skj-pop-tags_19820215.txt" file48="skj-pop-tags_19871215.txt" file49="skj-pop-tags_19880215.txt" file50="skj-pop-tags_19880315.txt" file51="skj-pop-tags_19880415.txt" file52="skj-pop-tags_19890215.txt" file53="skj-pop-tags_19890315.txt" file54="skj-pop-tags_19890415.txt" file55="skj-pop-tags_19890515.txt" file56="skj-pop-tags_19890615.txt" file57="skj-pop-tags_19900515.txt" file58="skj-pop-tags_19900915.txt" file59="skj-pop-tags_19910615.txt" file60="skj-pop-tags_19910715.txt" file61="skj-pop-tags_19920315.txt" file62="skj-pop-tags_19921015.txt" file63="skj-pop-tags_19930315.txt" file64="skj-pop-tags_19930415.txt" file65="skj-pop-tags_19930515.txt" file66="skj-pop-tags_19930815.txt" file67="skj-pop-tags_19930915.txt" file68="skj-pop-tags_19931015.txt" file69="skj-pop-tags_19931115.txt" file70="skj-pop-tags_19940415.txt" file71="skj-pop-tags_19940515.txt" file72="skj-pop-tags_19940615.txt" file73="skj-pop-tags_19940715.txt" file74="skj-pop-tags_19940815.txt" file75="skj-pop-tags_19940915.txt" file76="skj-pop-tags_19941015.txt" file77="skj-pop-tags_19950115.txt" file78="skj-pop-tags_19950215.txt" file79="skj-pop-tags_19950315.txt" file80="skj-pop-tags_19950415.txt" file81="skj-pop-tags_19950515.txt" file82="skj-pop-tags_19950615.txt" file83="skj-pop-tags_19950715.txt" file84="skj-pop-tags_19960215.txt" file85="skj-pop-tags_19960415.txt" file86="skj-pop-tags_19960515.txt" file87="skj-pop-tags_19970415.txt" file88="skj-pop-tags_19970515.txt" file89="skj-pop-tags_19970615.txt" file90="skj-pop-tags_19970715.txt" file91="skj-pop-tags_19990715.txt" file92="skj-pop-tags_19990915.txt" file93="skj-pop-tags_19991015.txt" file94="skj-pop-tags_20000515.txt" file95="skj-pop-tags_20000615.txt" file96="skj-pop-tags_20020615.txt" file97="skj-pop-tags_20040515.txt" file98="skj-pop-tags_20040615.txt" file99="skj-pop-tags_20060615.txt" file100="skj-pop-tags_20060815.txt" file101="skj-pop-tags_20060915.txt" file102="skj-pop-tags_20061015.txt" file103="skj-pop-tags_20061115.txt" file104="skj-pop-tags_20061215.txt" file105="skj-pop-tags_20070115.txt" file106="skj-pop-tags_20070215.txt" file107="skj-pop-tags_20070315.txt" file108="skj-pop-tags_20070415.txt" file109="skj-pop-tags_20070515.txt" file110="skj-pop-tags_20070615.txt" file111="skj-pop-tags_20070715.txt" file112="skj-pop-tags_20070815.txt" file113="skj-pop-tags_20070915.txt" file114="skj-pop-tags_20071015.txt" file115="skj-pop-tags_20071115.txt" file116="skj-pop-tags_20071215.txt" file117="skj-pop-tags_20080115.txt" file118="skj-pop-tags_20080215.txt" file119="skj-pop-tags_20080315.txt" file120="skj-pop-tags_20080415.txt" file121="skj-pop-tags_20080515.txt" file122="skj-pop-tags_20080615.txt" file123="skj-pop-tags_20080715.txt" file124="skj-pop-tags_20080815.txt" file125="skj-pop-tags_20080915.txt" file126="skj-pop-tags_20081015.txt" file127="skj-pop-tags_20081115.txt" file128="skj-pop-tags_20081215.txt" file129="skj-pop-tags_20090115.txt" file130="skj-pop-tags_20090215.txt" file131="skj-pop-tags_20090315.txt" file132="skj-pop-tags_20090415.txt" file133="skj-pop-tags_20090515.txt" file134="skj-pop-tags_20090615.txt" file135="skj-pop-tags_20090715.txt" file136="skj-pop-tags_20090815.txt" file137="skj-pop-tags_20090915.txt" file138="skj-pop-tags_20091015.txt" file139="skj-pop-tags_20091115.txt" file140="skj-pop-tags_20091215.txt" file141="skj-pop-tags_20100115.txt" file142="skj-pop-tags_20100215.txt" file143="skj-pop-tags_20100315.txt" file144="skj-pop-tags_20100415.txt" file145="skj-pop-tags_20100515.txt" file146="skj-pop-tags_20100615.txt" file147="skj-pop-tags_20100715.txt" file148="skj-pop-tags_20100815.txt" file149="skj-pop-tags_20100915.txt" file150="skj-pop-tags_20101015.txt" file151="skj-pop-tags_20101115.txt" file152="skj-pop-tags_20101215.txt" />
+</file_tag_data>
+<use_tag_masks flag="1"/>
+<strdir_tagmasks value="/home/lucasb/SEAPODYM/data/tagging/skj/ds2016/tagmasks/run-interim_2x30d_po/biomass5pcts/0.5knots/"/>
+
+<!--type of oxygen data: 0 (default) for time series, 1 for monthly, 2 for seasonal climatology-->
+
+<type_oxy value="0"/>
+
+<strfile_pp value="po_interim_historic_2x30d_pp_1979_2010.dym"/>
+<strfile_sst value="po_interim_historic_2x30d_sst_1979_2010.dym"/>
+<strfile_vld value="po_interim_historic_2x30d_zeu_1979_2010.dym"/>
+<strfile_u layer0="po_interim_historic_2x30d_u_L1_1979_2010.dym" layer1="po_interim_historic_2x30d_u_L2_1979_2010.dym" layer2="po_interim_historic_2x30d_u_L3_1979_2010.dym"/>
+<strfile_v layer0="po_interim_historic_2x30d_v_L1_1979_2010.dym" layer1="po_interim_historic_2x30d_v_L2_1979_2010.dym" layer2="po_interim_historic_2x30d_v_L3_1979_2010.dym"/>
+<strfile_t layer0="po_interim_historic_2x30d_temperature_L1_1979_2010.dym" layer1="po_interim_historic_2x30d_temperature_L2_1979_2010.dym" layer2="po_interim_historic_2x30d_temperature_L3_1979_2010.dym"/>
+<strfile_oxy layer0="po_interim_historic_2x30d_O2_L1_1979_2010.dym" layer1="po_interim_historic_2x30d_O2_L2_1979_2010.dym" layer2="po_interim_historic_2x30d_O2_L3_1979_2010.dym"/>
+
+<strfile_pp_clm value="po_interim_historic_2x30d_PP_1980.dym"/>
+<strfile_u_clm layer0="po_interim_historic_2x30d_u_L1_1980.dym" layer1="po_interim_historic_2x30d_u_L2_1980.dym" layer2="po_interim_historic_2x30d_u_L3_1980.dym"/>
+<strfile_v_clm layer0="po_interim_historic_2x30d_v_L1_1980.dym" layer1="po_interim_historic_2x30d_v_L2_1980.dym" layer2="po_interim_historic_2x30d_v_L3_1980.dym"/>
+<strfile_t_clm layer0="po_interim_historic_2x30d_temperature_L1_1980.dym" layer1="po_interim_historic_2x30d_temperature_L2_1980.dym" layer2="po_interim_historic_2x30d_temperature_L3_1980.dym"/>
+<strfile_oxy_clm layer0="po_interim_historic_2x30d_O2_L1_1980.dym" layer1="po_interim_historic_2x30d_O2_L2_1980.dym" layer2="po_interim_historic_2x30d_O2_L3_1980.dym"/>
+
+
+<!-- Output Directory -->
+<strdir_output value="output/"/>
+<write_all_cohorts_dym value="0"/>
+
+<!--====================== -->
+<!--FORAGE MODEL PARAMETERS-->
+<!--====================== -->
+<Tr_max value="450"/>
+<Tr_exp value="-0.1"/>
+<inv_lambda_max value="1825"/>
+<inv_lambda_curv value="-0.1"/>
+<E value="0.00316"/>
+<c_pp value="0.0948"/>
+
+<nb_forage value="6"/>
+
+<frg_name>epi meso mmeso bathy mbathy hmbathy</frg_name>
+
+<source_frg epi="0.17" meso="0.1" mmeso="0.22" bathy="0.18" mbathy="0.13" hmbathy="0.2"/>
+<day_layer epi="0" meso="1" mmeso="1" bathy="2" mbathy="2" hmbathy="2"/>
+<night_layer epi="0" meso="1" mmeso="0" bathy="2" mbathy="1" hmbathy="0"/>
+
+<tstep_forage value="30"/>
+<sigma_fcte value="10"/>
+<max_NPP value="200"/>
+<b_F_habitat value="20000"/>
+
+<!--=========================== -->
+<!--PREDATOR MODEL AGE STRUCTURE-->
+<!--=========================== -->
+<nb_species value="1"/>
+
+<sp_name>skj</sp_name>
+
+<nb_life_stages skj="3"/>
+<life_stage>
+        <skj>larvae juvenile adult</skj>
+</life_stage>
+
+<nb_cohort_life_stage>
+        <skj>1 2 34</skj>
+</nb_cohort_life_stage>
+
+<age_recruit skj="3"/>
+<age_mature skj="11"/>
+
+<sp_unit_cohort>
+    <skj> 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 30 600</skj>
+</sp_unit_cohort>
+<!--age = 0.5  1.5  2.5  3.5  4.5  5.5  6.5  7.5  8.5  9.5 10.5 11.5 12.5 13.5 14.5 15.5 16.5 17.5 18.5 19.5 20.5 21.5 22.5 23.5 24.5 25.5 26.5 27.5 28.5 29.5 30.5 31.5 32.5 33.5 34.5 35.5 36.5 37.5 38.5 39.5 40.5 41.5 42.5 43.5 44.5 45.5 46.5 47.5 -->
+<age_compute_habitat>
+  <skj>0  1  2  3  4  5  6  7  8  9 10 11 12 12 13 13 14 14 15 15 16 16 17 17 17
+       18 18 18 19 19 19 19 20 20 20 20 21 
+  </skj>
+</age_compute_habitat>
+
+<length>
+    <skj>3 4.5 7.52 11.65 16.91 21.83 26.43 30.72 34.73 38.49 41.99 45.27 48.33 51.19 53.86 56.36 58.70 60.88 62.92 64.83 66.61 68.27 69.83 71.28 72.64 73.91 75.10 76.21 77.25 78.22 79.12 79.97 80.76 81.50 82.19 82.83 87.96 </skj>
+<!--    <bounds>8.89 14.33 19.41 24.17 28.61 32.76 36.64 40.27 43.66 46.82 49.78 52.55 55.13 57.55 59.81 61.92 63.89 65.73 67.46 69.07 70.57 71.98 73.29 74.52 75.67 76.74 77.74 78.68 79.55 80.37 81.13 81.85 82.52 83.44 89.75 </bounds>-->
+</length>
+<weight>
+    <skj>0.0003 0.00109 0.00571 0.023 0.08  0.18  0.32  0.53  0.78  1.09  1.44  1.84 2.27  2.73  3.21  3.72  4.23  4.76  5.30  5.83  6.36  6.89  7.40  7.91 8.41  8.89  9.36  9.81 10.25 10.66 11.07 11.45 11.82 12.17 12.51 12.83 15.56  </skj>
+</weight>
+
+
+
+<!--======================== -->
+<!--PREDATOR MODEL PARAMETERS-->
+<!--======================== -->
+
+<!--=== 1. Available OPTIONS =====-->
+
+<seasonal_migrations skj="0"/>
+
+<vertical_movement skj="1"/>
+<!-- if optimal sst for larvae is assumed 
+     different from b_sst_spawning -->
+<uncouple_sst_larvae skj="1"/>
+
+<gaussian_thermal_function skj="1"/>
+
+<cannibalism_juv skj="0"/>
+
+<food_requirement_in_mortality skj="0"/>
+
+<!--= end of Available OPTIONS ==-->
+
+<!--===  2. Spawning habitat parameters  =======-->
+<a_sst_larvae skj="1.492626509052383">
+        <variable min="0.5" max="2.5" use="false"/>
+</a_sst_larvae>
+
+<b_sst_larvae skj="28.50028833256241"> 
+        <variable min="28.5" max="30.5" use="false"/>
+</b_sst_larvae>
+
+<a_sst_spawning skj="1.492626509052383">
+        <variable min="0.5" max="2.5" use="true"/>
+</a_sst_spawning>
+
+<b_sst_spawning skj="31.09996321152416">
+        <variable min="28.0" max="31.1" use="true"/>
+</b_sst_spawning>
+
+<alpha_hsp_prey skj="0.000569603615391101">
+        <variable min="0.0" max="1.0" use="false"/>
+</alpha_hsp_prey>
+
+<alpha_hsp_predator skj="2.30637293118113">
+        <variable min="0.0" max="2.5" use="false"/>
+</alpha_hsp_predator>
+
+<beta_hsp_predator skj="1.534250459427564">
+	<variable min="0.0" max="3.0" use="false"/>
+</beta_hsp_predator>
+<!--=== end of Spawning habitat parameters =====-->
+
+<!--===  3. Recruitment function   =============-->
+<nb_recruitment skj="0.03872619045606616">
+	<variable min="0.005" max="0.25" use="false"/>
+</nb_recruitment>
+
+<a_adults_spawning skj="0.1605221815904466">
+	<variable min="0.005" max="0.5" use="false"/>
+</a_adults_spawning>
+<!--=== end of Recruitment function  ===========-->
+
+<!--=== 4. Juvenile habitat parameters =========-->
+<hp_cannibalism skj="0.003505451231040554">
+	<variable min="0.002" max="0.3" use="false"/>
+</hp_cannibalism>
+<!--=== end of Juvenile habitat parameters =====-->
+
+<!--=== 5. Adult seasonal migrations ===========-->
+<spawning_season_peak skj="0.00323472646805989">
+        <variable min="0" max="40" use="false"/>
+</spawning_season_peak>
+
+<spawning_season_start skj="1.250054625443927">
+	<variable min="1.0" max="1.5" use="false"/>
+        </spawning_season_start>
+<!--=== end of Adult seasonal migrations =======-->
+
+<a_sst_habitat skj="3.000002576905231">
+        <variable min="3.0" max="6.5" use="true"/>
+        </a_sst_habitat>
+
+<!--=== 6. Adult feeding habitat ===============-->
+<b_sst_habitat skj="26.04960048814059">
+        <variable min="18.0" max="26.05" use="true"/>
+        </b_sst_habitat>
+
+<T_age_size_slope skj="2.83057972721579">
+        <variable min="0.5" max="3.0" use="true"/>
+</T_age_size_slope>
+
+<thermal_func>
+    <delta1 skj="0.04">
+        <variable min="0.0" max="0.05" use="false"/>
+    </delta1>
+
+    <delta2 skj="0.015">
+        <variable min="0.0" max="0.25" use="false"/>
+    </delta2>
+
+    <delta3 skj="0.25">
+        <variable min="0.0" max="10.0" use="false"/>
+    </delta3>
+</thermal_func>
+
+<a_oxy_habitat skj="0.0001000000005809866">
+        <variable min="0.0001" max="0.1" use="true"/>
+        </a_oxy_habitat>
+
+<b_oxy_habitat skj="3.764163854561827">
+        <variable min="1.0" max="4.0" use="true"/>
+        </b_oxy_habitat>
+
+<eF_habitat>
+    <epi skj="3.999907116416026">
+        <variable min="0.0" max="4.0" use="true"/>
+    </epi>
+
+    <meso skj="0.05">
+        <variable min="0.0" max="4.0" use="false"/>
+    </meso>
+
+    <mmeso skj="0.897307278945509">
+        <variable min="0.0" max="2.0" use="true"/>
+    </mmeso>
+
+    <bathy skj="0.0">
+        <variable min="0.0" max="1.0" use="false"/>
+    </bathy>
+
+    <mbathy skj="0.0">
+        <variable min="0.0" max="1.0" use="false"/>
+    </mbathy>
+
+    <hmbathy skj="3.999946416753347">
+        <variable min="0.0" max="4.0" use="true"/>
+    </hmbathy>
+</eF_habitat>
+<!--=== end of Adult feeding habitat ===========-->
+
+<!--=== 7. Movement parameters    ==============-->
+
+<sigma_species skj="4.805895">
+        <variable min="0.0001" max="5" use="true"/>
+        </sigma_species>
+
+<c_diff_fish skj="0.9294158175770607">
+        <variable min="0.1" max="1.0" use="true"/>
+        </c_diff_fish>
+
+<MSS_species skj="0.4893611123606742">
+        <variable min="0.0" max="5.0" use="true"/>
+        </MSS_species>
+	
+<MSS_size_slope skj="0.4278765609783742">
+        <variable min="0.4" max="1.15" use="true"/>
+        </MSS_size_slope>
+<!--=== end of Movement parameters    ==========-->
+
+<!--=== 8. Mortality parameters   ==============-->
+<Mp_mean_max skj="0.254521790396429624">
+        <variable min="0.0" max="0.5" use="false"/>
+</Mp_mean_max>
+
+<Mp_mean_exp skj="0.06595823270638622">
+        <variable min="0.05" max="0.55" use="false"/>
+</Mp_mean_exp>
+
+<Ms_mean_max skj="0.004521790396429624">
+        <variable min="0.0" max="0.01" use="false"/>
+</Ms_mean_max>
+
+<Ms_mean_slope skj="1.137368350830949">
+        <variable min="0.5" max="2.0" use="false"/>
+</Ms_mean_slope>
+
+<M_mean_range skj="2.30323700223245e-07">
+        <variable min="0.0" max="2.0" use="false"/>
+</M_mean_range>
+<!--=== end of Mortality parameters   ===========-->
+
+<!--=== 9. Food requirement index  ==============-->
+<!-- residual competition is the proportion of total F biomass 
+available for a given species in the presence of other species -->
+<residual_competition skj="0.025"/> 
+
+<!-- daily food ration -->
+<forage_ration skj="0.06"/>
+<!--=== end of Food requirement index  ==========-->
+
+
+<!--========================== -->
+<!--    FISHERY PARAMETERS     -->
+<!--========================== -->
+<eff_units_converter>1 1 1 1 1 1 1 1 1 1 1 1 1 1 1</eff_units_converter>
+<q_sp_fishery>
+
+	<variables use="true"/>
+
+	<P1 skj="0.001832370141253525">
+                <variable min="0.0" max="0.1" use="true" dyn="0.0"/>
+        </P1>
+
+        <P21 skj="0.0006686392664831038">
+                <variable min="0.0" max="0.05" use="false" dyn="0.0"/>
+        </P21>
+        
+        <P22 skj="0.001006072242453771">
+                <variable min="0.0" max="0.05" use="true" dyn="0.003"/>
+	</P22>
+
+        <P23 skj="0.001856932408223361">
+                <variable min="0.0" max="0.05" use="true" dyn="0.0"/>
+	</P23>
+
+	<P3 skj="0.0004184865115586356">
+                <variable min="0.0" max="0.05" use="true" dyn="0.0"/>
+        </P3>
+
+        <S4 skj="0.04177392903243383">
+                <variable min="0.0" max="0.2" use="false" dyn="0.0"/>
+        </S4>
+
+        <S5 skj="0.001802596303789339">
+                <variable min="0.0" max="0.2" use="false" dyn="0.0015"/>
+        </S5>
+
+        <S6 skj="0.004">
+                <variable min="0.0" max="0.1" use="false" dyn="0.0"/>
+        </S6>
+
+        <S7 skj="0.002071598118436371">
+                <variable min="0.0" max="0.05" use="false" dyn="0.0025"/>
+        </S7>
+
+        <L8 skj="0.0001424205818192381">
+                <variable min="0.0" max="0.01" use="false" dyn="0.0"/>
+        </L8>
+
+        <L9 skj="0.004">
+                <variable min="0.0" max="0.05" use="false" dyn="0.0"/>
+        </L9>
+
+        <S10 skj="0.00238844864774021">
+                <variable min="0.0" max="0.5" use="false" dyn="0.006"/>
+        </S10>
+
+        <S11 skj="0.01060339837354575">
+                <variable min="0.0" max="0.3" use="false" dyn="0.0"/>
+        </S11>
+
+        <S12 skj="0.0001570280274113101">
+                <variable min="0.0" max="0.2" use="false" dyn="0.0002"/>
+        </S12>
+
+        <S13 skj="0.001584022462382928">
+                <variable min="0.0" max="0.2" use="false" dyn="0.0035"/>
+        </S13>
+
+        <S14 skj="0.005357510497456042">
+                <variable min="0.0" max="0.2" use="false" dyn="0.0"/>
+        </S14>
+
+	<P15 skj="0.001757642365458623">
+                <variable min="0.0" max="0.05" use="true" dyn="0.0"/>
+        </P15>
+
+</q_sp_fishery>
+<!--selectivities as functions of length--> 
+<!--1 - logistic, 1 parameter: main parameter is used as a slope coefficient--> 
+<!--2 - sigmoid, 2 parameters: main parameter as slope coefficient + length_threshold--> 
+<!--3 - asymmetric Gaussian, 3 parameters: length_threshold as optimal value, main as sigma, right_asymptote as the function min on the right --> 
+<!--4 - two-modal Gaussian, 3 parameters: length_threshold as optimal value, main as sigma, right_asymptote as the second length --> 
+<s_sp_fishery>
+	<variables use="true"/>
+
+        <P1 skj="5.138119729769897">
+                <variable min="0.05" max="20.0" use="true"/>
+                <function_type value="3"/>
+                <length_threshold skj="48.75922407055393" use="true"/>
+                <right_asymptote skj="0.7051596508008535" min="0.0" max="1" use="true"/>
+        </P1>                                                                                                                                                            
+        <P21 skj="4.0">
+                <variable min="0.5" max="26.0" use="false"/>
+                <function_type value="3"/>
+                <length_threshold skj="37.0" use="false"/>
+                <right_asymptote skj="1.0" min="0.1" max="1" use="false"/>
+        </P21>                                                                                                                                                            
+        <P22 skj="10.0">
+                <variable min="1.0" max="26.0" use="false"/>
+                <function_type value="3"/>
+                <length_threshold skj="54.0" use="false"/>
+                <right_asymptote skj="0.9999999766914867" min="0.0" max="1" use="true"/>
+        </P22>                                                                                                                                                            
+        <P23 skj="12.79918705574716">
+                <variable min="1.0" max="26.0" use="true"/>
+                <function_type value="3"/>
+                <length_threshold skj="59.6204358140092" use="true"/>
+                <right_asymptote skj="0.9999999699539683" min="0.0" max="1" use="true"/>
+        </P23>                                                                                                                                                            
+        <P3 skj="9.5">
+                <variable min="0.05" max="20.0" use="false"/>
+                <function_type value="3"/>
+                <length_threshold skj="48.0" use="false"/>
+                <right_asymptote skj="0.05" min="0.0" max="1" use="false"/>
+        </P3>                                                                                                                                                            
+        <S4 skj="6">
+                <variable min="0.05" max="20.0" use="false"/>
+                <function_type value="3"/>
+                <length_threshold skj="52.0" use="false"/>
+                <right_asymptote skj="0" min="0.0" max="1" use="false"/>
+        </S4>
+
+        <S5 skj="12.6936723589497">
+                <variable min="0.05" max="20.0" use="true"/>
+                <function_type value="3"/>
+                <length_threshold skj="52.66207651171928" use="true"/>
+                <right_asymptote skj="0.732822821995287" min="0.0" max="1" use="true"/>
+        </S5>
+
+        <S6 skj="14.0">
+                <variable min="0.05" max="20.0" use="false"/>
+                <function_type value="3"/>
+                <length_threshold skj="64.0" use="false"/>
+                <right_asymptote skj="0.01" min="0.0" max="1" use="false"/>
+	</S6>
+
+        <S7 skj="16.26887151879928">
+                <variable min="0.05" max="24.0" use="true"/>
+                <function_type value="3"/>
+                <length_threshold skj="62.45871392820482" use="true"/>
+                <right_asymptote skj="0.9999993114821045" min="0.0" max="1" use="true"/>
+        </S7>   
+
+        <L8 skj="0.5">
+                <variable min="0.0" max="1" use="false"/>
+                <function_type value="2"/>
+                <length_threshold skj="67.0" use="false"/>
+                <right_asymptote skj="0.85" min="0" max="1" use="false"/>
+        </L8>                                                                                                                                                            
+        <L9 skj="0.3">
+                <variable min="0" max="1" use="false"/>
+                <function_type value="2"/>
+                <length_threshold skj="60" use="false"/>
+                <right_asymptote skj="0.8" min="0" max="1" use="false"/>
+        </L9>   
+
+        <S10 skj="11.31333403595283">
+                <variable min="0.5" max="26" use="true"/>
+                <function_type value="3"/>
+                <length_threshold skj="47.65341035786737" use="true"/>
+                <right_asymptote skj="0.9999999990945595" min="0.1" max="1" use="true"/>
+        </S10>                                                                                                                                                            
+
+        <S11 skj="19.35942376923553">
+                <variable min="0.5" max="30" use="true"/>
+                <function_type value="3"/>
+                <length_threshold skj="68.59420043098373" use="true"/>
+                <right_asymptote skj="0.9999998444525892" min="0.1" max="1" use="true"/>
+        </S11>                                                                                                                                                            
+
+        <S12 skj="14.97999590862408">
+                <variable min="0.5" max="26.0" use="true"/>
+                <function_type value="3"/>
+                <length_threshold skj="65.2624359755044" use="true"/>
+                <right_asymptote skj="0.9999988045096152" min="0.1" max="1" use="true"/>
+        </S12> 
+
+        <S13 skj="22.81650942479624">
+                <variable min="0.0" max="28" use="true"/>
+                <function_type value="3"/>
+                <length_threshold skj="71.38834887459542" use="true"/>
+                <right_asymptote skj="0.9999949108516821" min="0.1" max="1" use="true"/>
+        </S13>
+
+        <S14 skj="21.55910703359929">
+                <variable min="0.0" max="28" use="true"/>
+                <function_type value="3"/>
+                <length_threshold skj="71.60314451069418" use="true"/>
+                <right_asymptote skj="0.9999815195952257" min="0.0" max="1" use="true"/>
+        </S14>                                                                                                                                                                
+        <P15 skj="15.0">
+                <variable min="0.05" max="26.0" use="false"/>
+                <function_type value="3"/>
+                <length_threshold skj="65.0" use="false"/>
+                <right_asymptote skj="0.5" min="0.1" max="1" use="false"/>
+        </P15>                                                                                                                                                            
+ </s_sp_fishery>
+
+
+
+<!--================ -->
+<!--AGGREGATION ZONES-->
+<!--================ -->
+<!--Starting 2014 5 regions in WCPO with two given by composition of several rectangles, 
+     which is not implemented, so will declare simple, but modify MFCL file to provide
+     R2 = R2+R5 -->
+<nb_region value="5"/>
+<use_mask_catch value="1"/>
+<area0 area_id="1" lgmin="120" lgmax="210" ltmin="20" ltmax="50"/>
+<area1 area_id="2" lgmin="140" lgmax="170" ltmin="-20" ltmax="20"/>
+<area2 area_id="3" lgmin="170" lgmax="210" ltmin="-20" ltmax="20"/>
+<area3 area_id="4" lgmin="110" lgmax="140" ltmin="-20" ltmax="20"/>
+<area4 area_id="5" lgmin="210" lgmax="290" ltmin="-20" ltmax="30"/>
+
+<nb_region_sp_B skj="5"/>
+<area_sp_B> 
+	<skj>1 2 3 4 5</skj>
+</area_sp_B>
+<!--
+<nb_region value="7"/>
+<use_mask_catch value="0"/>
+<area0 area_id="1" lgmin="120.5" lgmax="140.5" ltmin="25.5" ltmax="35.5"/>
+<area1 area_id="2" lgmin="140.5" lgmax="165.5" ltmin="25.5" ltmax="45.5"/>
+<area2 area_id="3" lgmin="165.5" lgmax="210.5" ltmin="25.5" ltmax="45.5"/>
+<area3 area_id="4" lgmin="115.5" lgmax="210.5" ltmin="15.5" ltmax="25.5"/>
+<area4 area_id="5" lgmin="115.5" lgmax="165.5" ltmin="-20.5" ltmax="15.5"/>
+<area5 area_id="6" lgmin="165.5" lgmax="210.5" ltmin="-20.5" ltmax="15.5"/>
+<area6 area_id="7" lgmin="210.5" lgmax="280.5" ltmin="-20.5" ltmax="25.5"/>
+
+<nb_region_sp_B skj="7"/>
+<area_sp_B>
+        <skj>1 2 3 4 5 6 7 8 9</skj>
+</area_sp_B>
+-->
+<!--
+<nb_EEZ value="0"/>
+<str_file_maskEEZ value="mask_omega_PNG_AW.txt"/>
+<EEZ>
+        <eez0 name="NORTH" id="-1"/>
+        <eez1 name="Bismarck Sea" id="-2"/>
+        <eez2 name="EAST" id="-3"/>
+        <eez3 name="Solomon Sea" id="-4"/>
+</EEZ>
+-->
+<nb_EEZ value="0"/>
+<str_file_maskEEZ value="po_interim_1deg_mask_eez_wcpfc.txt"/>
+<EEZ>
+        <eez0 name="FSM" id="-13"/>
+        <eez1 name="Kiribati" id="-14"/>
+        <eez2 name="Marshall Islands" id="-17"/>
+        <eez3 name="Nauru" id="-21"/>
+        <eez4 name="Palau" id="-27"/>
+        <eez5 name="PNG" id="-25"/>
+        <eez6 name="Solomon Islands" id="-29"/>
+        <eez7 name="Tuvalu" id="-32"/>
+        <eez8 name="Cook Island" id="-11"/>
+</EEZ>
+
+<!--
+<nb_EEZ value="0"/>
+<str_file_maskEEZ value="mask_eez_mariana_islands.txt"/>
+<EEZ>
+        <eez0 name="Northern Mariana Islands" id="-2"/>
+        </EEZ>
+-->
+
+<likelihood_parameters>
+	<variables use="false"/>
+	<P1 skj="0.01"/>
+        <P21 skj="0.007"/>
+        <P22 skj="0.007"/>
+        <P23 skj="0.007"/>
+        <P3 skj="0.0"/>
+        <S4 skj="0.5"/>
+	<S5 skj="0.01"/>
+	<S6 skj="0.01"/>
+	<S7 skj="0.01"/>
+	<L8 skj="0.5"/>
+	<L9 skj="0.01"/>
+	<S10 skj="0.25"/>
+	<S11 skj="0.01"/>
+	<S12 skj="0.01"/>
+	<S13 skj="0.01"/>
+	<S14 skj="0.01"/>
+	<P15 skj="0.01"/>
+</likelihood_parameters>
+
+<prob_zero>
+	<P1 skj="0.0"/>
+        <P21 skj="0.0"/>
+        <P22 skj="0.0"/>
+        <P23 skj="0.0"/>
+        <P3 skj="0.0"/>
+        <S4 skj="0.0"/>
+	<S5 skj="0.0"/>
+	<S6 skj="0.0"/>
+	<S7 skj="0.0"/>
+	<L8 skj="0.0"/>
+	<L9 skj="0.0"/>
+	<S10 skj="0.0"/>
+	<S11 skj="0.0"/>
+	<S12 skj="0.0"/>
+	<S13 skj="0.0"/>
+	<S14 skj="0.0"/>
+        <P15 skj="0.0"/>
+</prob_zero>
+
+</par>