Working version of runmodelbyLME.R & dbpm_model_functions.r & PLottin…

…g_functions_DBPM.r - test LME 1 and gravity option 2 (issues from option 1 solved)
Benthic-Pelagic-Size-Spectrum-Model · Aug 24, 2023 · a0fbc79 · a0fbc79
1 parent fc22b58
commit a0fbc79
Show file tree

Hide file tree

Showing 13 changed files with 24 additions and 45 deletions.
diff --git a/Output/Effort_LME1_check.pdf b/Output/Effort_LME1_check.pdf
diff --git a/Output/TEST1_rungridbyLME.pdf → Output/Gravity0_biomass_LME1.pdf b/Output/TEST1_rungridbyLME.pdf → Output/Gravity0_biomass_LME1.pdf
diff --git a/Output/TEST1_catches_rungridbyLME.pdf → Output/Gravity0_catches_LME1.pdf b/Output/TEST1_catches_rungridbyLME.pdf → Output/Gravity0_catches_LME1.pdf
diff --git a/Output/Gravity1_BiomassU_LME_1_fishing.pdf b/Output/Gravity1_BiomassU_LME_1_fishing.pdf
diff --git a/Output/Gravity1_BiomassV_LME_1_fishing.pdf b/Output/Gravity1_BiomassV_LME_1_fishing.pdf
diff --git a/Output/Gravity1_CatchU_LME_1_fishing.pdf b/Output/Gravity1_CatchU_LME_1_fishing.pdf
diff --git a/Output/Gravity1_CatchV_LME_1_fishing.pdf b/Output/Gravity1_CatchV_LME_1_fishing.pdf
diff --git a/Output/Gravity2_BiomassU_LME_1_fishing.pdf b/Output/Gravity2_BiomassU_LME_1_fishing.pdf
diff --git a/Output/Gravity2_BiomassV_LME_1_fishing.pdf b/Output/Gravity2_BiomassV_LME_1_fishing.pdf
diff --git a/Output/Gravity2_CatchU_LME_1_fishing.pdf b/Output/Gravity2_CatchU_LME_1_fishing.pdf
diff --git a/Output/Gravity2_CatchV_LME_1_fishing.pdf b/Output/Gravity2_CatchV_LME_1_fishing.pdf
diff --git a/Plotting_functions_DBPM.R b/Plotting_functions_DBPM.R
@@ -299,32 +299,32 @@ plotgridbyLME<-function(LMEnumber = 1){
 
   # print also in local for fast checking 
   name = ifelse(f.effort == FALSE, "_no_fishing.pdf", "_fishing.pdf")
-  ## adding gravity option to plot name for testing 
-  pdf(paste0("Output/Gravity1_LME_",LMEnumber, name), height = 8, width = 6, onefile = T)
-  print(p1)
-  print(p2)
-  print(p3)
-  print(p4)
-  print(p5)
-  print(p6)
-  print(p7)
-  print(p8)
-  dev.off()
+  # ## adding gravity option to plot name for testing 
+  # pdf(paste0("Output/Gravity2_LME_",LMEnumber, name), height = 8, width = 6, onefile = T)
+  # print(p1)
+  # print(p2)
+  # print(p3)
+  # print(p4)
+  # print(p5)
+  # print(p6)
+  # print(p7)
+  # print(p8)
+  # dev.off()
 
   # increase the size of some plots as are not visible 
-  pdf(paste0("Output/Gravity1_BiomassU_LME_",LMEnumber, name))
+  pdf(paste0("Output/Gravity2_BiomassU_LME_",LMEnumber, name))
   p1
   dev.off()
 
-  pdf(paste0("Output/Gravity1_BiomassV_LME_",LMEnumber, name))
+  pdf(paste0("Output/Gravity2_BiomassV_LME_",LMEnumber, name))
   p2
   dev.off()
 
-  pdf(paste0("Output/Gravity1_CatchU_LME_",LMEnumber, name))
+  pdf(paste0("Output/Gravity2_CatchU_LME_",LMEnumber, name))
   p6
   dev.off()
 
-  pdf(paste0("Output/Gravity1_CatchV_LME_",LMEnumber, name))
+  pdf(paste0("Output/Gravity2_CatchV_LME_",LMEnumber, name))
   p7
   dev.off()
 

diff --git a/dbpm_model_functions.R b/dbpm_model_functions.R
@@ -282,33 +282,13 @@ gridded_sizemodel<-function(params,ERSEM.det.input=F,U_mat,V_mat,W_mat,temp.effe
 
         # get proportion of total fishable biomass for each grid cell
         #output rates of fisheries catches per yr at size
-
-        # CN check
-        # class(prop.b) # vector - grid cell
         prop.b<-(apply(U[,Fref.u:Nx,i]*10^x[Fref.u:Nx]*dx,1,sum) + apply(V[,Fref.v:Nx,i]*10^x[Fref.v:Nx]*dx,1,sum))/sum(apply(U[,Fref.u:Nx,i]*10^x[Fref.u:Nx]*dx,1,sum) + apply(V[,Fref.v:Nx,i]*10^x[Fref.v:Nx]*dx,1,sum))
         #check sums to 1
         #sum(prop.b)
 
         # redistribute total effort across gridcells according to proportion of biomass in that gridcell
-
-        # # CN check
-        # dim(effort) # grid cell X time
-        # effort[,3]
-
-        effort[,i+1] = gravitymodel(effort[,i+1], prop.b,depth = depth,iter=10000)
-
-        # sum(effort[,1])
-        # sum(effort[,2]) # OK same starting values so should be the same after spreading
-        #
-        # # CN check
-        # effort[,1]
-        # effort[,2]
-        # effort[,3]
-        #
-        # # check different gravity models:
-        # # option 2 # cannot check like this - need to check in terms of catches
-        # effort[,2]
-        # sum(effort[,2])
+        effort[,i+1] = gravitymodel(effort[,i+1], prop.b,depth = depth,iter=1)
+        # for option 2 iter = 1
 
       } # end gravity model 
 
@@ -558,20 +538,19 @@ gridded_sizemodel<-function(params,ERSEM.det.input=F,U_mat,V_mat,W_mat,temp.effe
 
       }		
 
-      #output rates of fisheries catches per yr at size
-
-      ### WARNING adding j to the below 
+      # update Fmort 
       if (i < Neq){
-      Fvec.u[j,Fref.u:Nx,i+1] = Fmort.u*effort[j,i+1]
-      Fvec.v[j,Fref.v:Nx,i+1] = Fmort.v*effort[j,i+1]
+        Fvec.u[j,Fref.u:Nx,i+1] = Fmort.u*effort[j,i+1]
+        Fvec.v[j,Fref.v:Nx,i+1] = Fmort.v*effort[j,i+1]
       }
 
-      #### WARNING adding j to Y.u and Y.v 
+      #output rates of fisheries catches per yr at size
       Y.u[j,Fref.u:Nx,i+1]<-Fvec.u[j,Fref.u:Nx,i+1]*U[j,Fref.u:Nx,i+1]*10^x[Fref.u:Nx]
       #output rates of fisheries catches per yr at size
       Y.v[j,Fref.v:Nx,i+1]<-Fvec.v[j,Fref.v:Nx,i+1]*V[j,Fref.v:Nx,i+1]*10^x[Fref.v:Nx]
 
       # ### ORIGINAL
+      #output rates of fisheries catches per yr at size
       # Y.u[,Fref.u:Nx,i+1]<-Fvec.u[j,Fref.u:Nx,i]*U[j,Fref.u:Nx,i+1]*10^x[Fref.u:Nx]
       # #output rates of fisheries catches per yr at size
       # Y.v[,Fref.v:Nx,i+1]<-Fvec.v[j,Fref.v:Nx,i]*V[j,Fref.v:Nx,i+1]*10^x[Fref.v:Nx]
@@ -1484,9 +1463,9 @@ gravitymodel<-function(effort=effort[,3000],prop.b, depth, iter){
   for(j in 1:iter) {
 
     a<-eff
-    suit = prop.b*(1-d/max(d))*(1-a/0.001)
+    # suit = prop.b*(1-d/max(d))*(1-a/0.001)
     # # CN option 2 from below 
-    # suit = prop.b*(1-d/max(d))
+    suit = prop.b*(1-d/max(d))
     # # CN option 3 from below 
     # suit = prop.b
     # rescale: