From 253f0e2cd6096f82eaa966000062dc65d80e5893 Mon Sep 17 00:00:00 2001
From: KyleKlenk <kyle.c.klenk@gmail.com>
Date: Thu, 5 May 2022 13:27:39 -0600
Subject: [PATCH] Writing without Segfaulting

---
 build/source/actors/FileAccess.h              |  4 +-
 build/source/actors/FileAccessActor.h         | 70 ++++++++++++++----
 build/source/actors/HRUActor.h                |  2 -
 build/source/actors/OutputManager.h           | 65 ++++++++++++++---
 .../file_access_actor/cppwrap_fileAccess.f90  | 68 +++++++++--------
 build/source/netcdf/writeOutput.f90           | 73 ++++++++++---------
 6 files changed, 189 insertions(+), 93 deletions(-)

diff --git a/build/source/actors/FileAccess.h b/build/source/actors/FileAccess.h
index 3650ea5..063dc3d 100644
--- a/build/source/actors/FileAccess.h
+++ b/build/source/actors/FileAccess.h
@@ -5,6 +5,7 @@
 #include "../interface/file_access_actor/fileAccess_subroutine_wrappers.h"
 #include "caf/all.hpp"
 #include "messageAtoms.h"
+#include "OutputManager.h"
 #include <vector>
 #include <chrono>
 #include "global.h"
@@ -48,7 +49,8 @@ struct file_access_state {
 
     void *handle_forcFileInfo = new_handle_file_info(); // Handle for the forcing file information
     void *handle_ncid = new_handle_var_i();               // output file ids
-    bool outputFileExists = false;
+    OutputManager *output_manager;
+    int num_vectors_in_output_manager = 5;
     int num_steps;
     int outputStrucSize;
     int stepsInCurrentFile;
diff --git a/build/source/actors/FileAccessActor.h b/build/source/actors/FileAccessActor.h
index c3d1d59..6f30888 100644
--- a/build/source/actors/FileAccessActor.h
+++ b/build/source/actors/FileAccessActor.h
@@ -5,7 +5,7 @@
 
 using namespace caf;
 void initalizeFileAccessActor(stateful_actor<file_access_state>* self);
-int writeOutput(stateful_actor<file_access_state>* self, int indxGRU, int indxHRU, int numStepsToWrite);
+int writeOutput(stateful_actor<file_access_state>* self, int indxGRU, int indxHRU, int numStepsToWrite, int returnMessage, caf::actor actorRef);
 int readForcing(stateful_actor<file_access_state>* self, int currentFile);
 
 behavior file_access_actor(stateful_actor<file_access_state>* self, int startGRU, int numGRU, 
@@ -102,13 +102,12 @@ behavior file_access_actor(stateful_actor<file_access_state>* self, int startGRU
         [=](write_output, int indxGRU, int indxHRU, int numStepsToWrite,
             caf::actor refToRespondTo) {
             int err;
+            int returnMessage = 9999;
             
-            err = writeOutput(self, indxGRU, indxHRU, numStepsToWrite);
+            err = writeOutput(self, indxGRU, indxHRU, numStepsToWrite, returnMessage, refToRespondTo);
             if (err != 0) {
                 aout(self) << "FILE_ACCESS_ACTOR - ERROR Writing Output \n";
-            } else {
-                self->send(refToRespondTo, done_write_v);
-            }
+            } 
 
 
         },
@@ -117,17 +116,13 @@ behavior file_access_actor(stateful_actor<file_access_state>* self, int startGRU
             caf::actor refToRespondTo) {
             int err;
 
-            err = writeOutput(self, indxGRU, indxHRU, numStepsToWrite);
-            if (err != 0)
-                aout(self) << "FILE_ACCESS_ACTOR - ERROR Writing Output \n";
-
             err = readForcing(self, currentFile);
             if (err != 0)
                 aout(self) << "\nERROR: FILE_ACCESS_ACTOR - READING_FORCING FAILED\n";
 
-            // Respond to HRU
-            self->send(refToRespondTo, run_hru_v, 
-                self->state.forcFileList[currentFile - 1].getNumSteps());
+            err = writeOutput(self, indxGRU, indxHRU, numStepsToWrite, currentFile, refToRespondTo);
+            if (err != 0)
+                aout(self) << "FILE_ACCESS_ACTOR - ERROR Writing Output \n";
         },
 
         [=](deallocate_structures) {
@@ -205,8 +200,10 @@ void initalizeFileAccessActor(stateful_actor<file_access_state>* self) {
     aout(self) << "Initalizing Output Structure" << std::endl;
     Init_OutputStruct(self->state.handle_forcFileInfo, &self->state.outputStrucSize, 
         &self->state.numGRU, &self->state.err);
-    
 
+    // Initalize the output manager  
+    self->state.output_manager = new OutputManager(self->state.num_vectors_in_output_manager, self->state.numGRU);
+    
     self->send(self->state.parent, done_file_access_actor_init_v);
     // initalize the forcingFile array
     self->state.filesLoaded = 0;
@@ -216,12 +213,53 @@ void initalizeFileAccessActor(stateful_actor<file_access_state>* self) {
 }
 
 int writeOutput(stateful_actor<file_access_state>* self, int indxGRU, int indxHRU, 
-    int numStepsToWrite) {
+    int numStepsToWrite, int returnMessage, caf::actor actorRef) {
     self->state.writeStart = std::chrono::high_resolution_clock::now();
+    if (debug) {
+        aout(self) << "Recieved Write Request From GRU: " << indxGRU << "\n";
+    }
+
+    
     int err = 0;
+    int listIndex = self->state.output_manager->addActor(actorRef, indxGRU, returnMessage);
+    if (self->state.output_manager->isFull(listIndex)) {
+        if (debug) {
+            aout(self) << "List with Index " << listIndex << " is full and ready to write\n";
+            aout(self) << "Minimum GRU Index = " << self->state.output_manager->getMinIndex(listIndex) << "\n";
+            aout(self) << "Maximum GRU Index = " << self->state.output_manager->getMaxIndex(listIndex) << "\n";
+        }
+
+        int minGRU = self->state.output_manager->getMinIndex(listIndex);
+        int maxGRU = self->state.output_manager->getMaxIndex(listIndex);
+        FileAccessActor_WriteOutput(self->state.handle_ncid,
+            &numStepsToWrite, &minGRU, 
+            &maxGRU, &err);
+        
+        // Pop The actors and send them the correct continue message
+        while(!self->state.output_manager->isEmpty(listIndex)) {
+            std::tuple<caf::actor, int> actor = self->state.output_manager->popActor(listIndex);
+            if (get<1>(actor) == 9999) {
+                
+                self->send(get<0>(actor), done_write_v);
+
+            }  else {
+                self->send(get<0>(actor), run_hru_v, 
+                    self->state.forcFileList[get<1>(actor) - 1].getNumSteps());
+            }
+        }
+        
+
+    } else {
+        if (debug) {
+            aout(self) << "List with Index " << listIndex << " is not full yet waiting to write\n";
+            aout(self) << "Size of list is " << self->state.output_manager->getSize(listIndex) << "\n";
+        }
+    }
     
-    FileAccessActor_WriteOutput(self->state.handle_ncid,
-        &numStepsToWrite, &indxGRU, &indxHRU, &err);
+   
+
+
+
     self->state.writeEnd = std::chrono::high_resolution_clock::now();
     self->state.writeDuration += calculateTime(self->state.writeStart, self->state.writeEnd);
 
diff --git a/build/source/actors/HRUActor.h b/build/source/actors/HRUActor.h
index 81b085a..a0dff80 100644
--- a/build/source/actors/HRUActor.h
+++ b/build/source/actors/HRUActor.h
@@ -455,8 +455,6 @@ bool check_HRU(stateful_actor<hru_state>* self, int err) {
         return false;
 
     } else {
-        if (debug)
-            aout(self) << "Continuing\n";
         return true;
     }
 }
diff --git a/build/source/actors/OutputManager.h b/build/source/actors/OutputManager.h
index af99740..aa3fa4c 100644
--- a/build/source/actors/OutputManager.h
+++ b/build/source/actors/OutputManager.h
@@ -3,12 +3,17 @@
 
 #include "caf/all.hpp"
 #include <vector>
-
+/**
+ * @brief Basic Container class to hold actor references. This has a size component for checking when it is full.
+ * 
+ */
 class ActorRefList {
     private:
         int currentSize;
         int maxSize;
-        std::vector<caf::actor> list;
+        int minIndex = -1; // minimum index of the actor being stored on this list
+        int maxIndex = 0; // maximum index of the actor being stored on this list
+        std::vector<std::tuple<caf::actor, int>> list;
     
     public:
         // Constructor
@@ -19,6 +24,14 @@ class ActorRefList {
 
         // Deconstructor
         ~ActorRefList(){};
+
+        int getMaxIndex() {
+            return this->maxIndex;
+        }
+
+        int getMinIndex() {
+            return this->minIndex;
+        }
         
         int getCurrentSize() {
             return this->currentSize;
@@ -28,17 +41,24 @@ class ActorRefList {
             return list.size() == this->maxSize;
         }
 
-        void addActor(caf::actor actor) {
+        void addActor(caf::actor actor, int index, int returnMessage) {
             if (this->isFull()) {
                 throw "List is full, cannot add actor to this list";
             }
+            if (index > this->maxIndex) {
+                this->maxIndex = index;
+            } 
+            if (index < this->minIndex || this->minIndex < 0) {
+                this->minIndex = index;
+            }
+
             this->currentSize++;
-            list.push_back(actor);
+            list.push_back(std::make_tuple(actor, returnMessage));
         }
 
-        caf::actor popActor() {
+        std::tuple<caf::actor,int> popActor() {
             if (list.empty()) {
-                throw "List is empty, nothing to pop.";
+                throw "List is empty, nothing to pop";
             }
             auto actor = list.back();
             list.pop_back();
@@ -46,6 +66,7 @@ class ActorRefList {
             return actor;
         }
 
+
         bool isEmpty() {
             return list.empty();
         }
@@ -55,6 +76,11 @@ class ActorRefList {
 
 };
 
+
+/**
+ * @brief Class that manages which structure actors are held on
+ * 
+ */
 class OutputManager {
     private:
 
@@ -83,17 +109,25 @@ class OutputManager {
         // Deconstructor
         ~OutputManager(){};
 
-        void addActor(caf::actor actor, int index) {
+        /**
+         * @brief Adds an actor to its respective list
+         * 
+         * @param actor Actor reference
+         * @param index Actor Index
+         * @return int The list index that actor is added to.
+         */
+        int addActor(caf::actor actor, int index, int returnMessage) {
             // Index has to be subtracted by 1 because Fortran array starts at 1
             int listIndex = (index - 1) / this->avgSizeOfActorList;
             if (listIndex > this->numVectors - 1) {
                 listIndex =  this->numVectors - 1;
             }
 
-            this->list[listIndex]->addActor(actor);
+            this->list[listIndex]->addActor(actor, index, returnMessage);
+            return listIndex;
         }
 
-        caf::actor popActor(int index) {
+        std::tuple<caf::actor,int> popActor(int index) {
             if (index > this->numVectors - 1 || index < 0) {
                 throw "List Index Out Of Range";
             } else if (this->list[index]->isEmpty()) {
@@ -103,6 +137,7 @@ class OutputManager {
             return this->list[index]->popActor();
 
         }
+        
 
         bool isFull(int listIndex) {
             if (listIndex > this->numVectors - 1) {
@@ -111,6 +146,10 @@ class OutputManager {
             return this->list[listIndex]->isFull();
         }
 
+        bool isEmpty(int listIndex) {
+            return this->list[listIndex]->isEmpty();
+        }
+
         int getSize(int listIndex) {
             if (listIndex > this->numVectors - 1) {
                 throw "List Index Out Of Range";
@@ -118,6 +157,14 @@ class OutputManager {
             return this->list[listIndex]->getCurrentSize();
         }
 
+        int getMinIndex(int listIndex) {
+            return this->list[listIndex]->getMinIndex();
+        }
+
+        int getMaxIndex(int listIndex) {
+            return this->list[listIndex]->getMaxIndex();
+        }
+
 };
 
 #endif 
\ No newline at end of file
diff --git a/build/source/interface/file_access_actor/cppwrap_fileAccess.f90 b/build/source/interface/file_access_actor/cppwrap_fileAccess.f90
index 4bcd544..db39100 100644
--- a/build/source/interface/file_access_actor/cppwrap_fileAccess.f90
+++ b/build/source/interface/file_access_actor/cppwrap_fileAccess.f90
@@ -279,10 +279,10 @@ subroutine Write_HRU_Param(&
 end subroutine
 
 subroutine FileAccessActor_WriteOutput(&
-                                handle_ncid,      & ! ncid of the output file
-                                nSteps,           & ! number of steps to write
-                                indxGRU,          & ! index of GRU we are currently writing for
-                                indxHRU,          & ! index of HRU we are currently writing for
+                                handle_ncid,     & ! ncid of the output file
+                                nSteps,          & ! number of steps to write
+                                minGRU,          & ! index of GRU we are currently writing for
+                                maxGRU,          & ! index of HRU we are currently writing for
                                 err) bind(C, name="FileAccessActor_WriteOutput")
   USE def_output_module,only:def_output                       ! module to define model output
   USE globalData,only:gru_struc
@@ -302,8 +302,8 @@ subroutine FileAccessActor_WriteOutput(&
   ! dummy variables
   type(c_ptr),intent(in), value        :: handle_ncid       ! ncid of the output file
   integer(c_int),intent(in)            :: nSteps            ! number of steps to write
-  integer(c_int),intent(in)            :: indxGRU           ! index of GRU we are currently writing for
-  integer(c_int),intent(in)            :: indxHRU           ! index of HRU we are currently writing for
+  integer(c_int),intent(in)            :: minGRU           ! index of GRU we are currently writing for
+  integer(c_int),intent(in)            :: maxGRU           ! index of HRU we are currently writing for
   integer(c_int),intent(inout)         :: err               ! Error code
   
   ! local variables 
@@ -311,56 +311,66 @@ subroutine FileAccessActor_WriteOutput(&
   character(LEN=256)                   :: message
   character(LEN=256)                   :: cmessage
   integer(i4b)                         :: iStruct
+  integer(i4b)                         :: numGRU
+  integer(i4b)                         :: iGRU
   integer(i4b)                         :: iStep
   integer(i4b)                         :: iFreq
+  integer(i4b)                         :: indxHRU=1
   integer(i4b), dimension(maxVarFreq)  :: outputTimestepUpdate
 
   call c_f_pointer(handle_ncid, ncid)
   ! ****************************************************************************
   ! *** write data
   ! ****************************************************************************
-  do iStep=1, nSteps
-    call writeBasin(ncid,indxGRU,outputTimeStep(indxGRU)%dat(:),iStep,bvar_meta, &
-              outputStructure(1)%bvarStat(1)%gru(indxGRU)%hru(indxHRU)%var, &
-              outputStructure(1)%bvarStruct(1)%gru(indxGRU)%hru(indxHRU)%var, bvarChild_map, err, cmessage)
+  do iGRU=minGRU, maxGRU
+    do iStep=1, nSteps
+      call writeBasin(ncid,iGRU,outputTimeStep(iGRU)%dat(:),iStep,bvar_meta, &
+              outputStructure(1)%bvarStat(1)%gru(iGRU)%hru(indxHRU)%var, &
+              outputStructure(1)%bvarStruct(1)%gru(iGRU)%hru(indxHRU)%var, bvarChild_map, err, cmessage)
   
 
-    ! reset outputTimeStep
-    ! get the number of HRUs in the run domain
-    ! nHRUrun = sum(gru_struc%hruCount)
-    ! write time information
-    call writeTime(ncid,outputTimeStep(indxGRU)%dat(:),iStep,time_meta, &
-              outputStructure(1)%timeStruct(1)%gru(indxGRU)%hru(indxHRU)%var,err,cmessage)
-  end do ! istep
-  
+      ! reset outputTimeStep
+      ! get the number of HRUs in the run domain
+      ! nHRUrun = sum(gru_struc%hruCount)
+      ! write time information
+      call writeTime(ncid,outputTimeStep(iGRU)%dat(:),iStep,time_meta, &
+              outputStructure(1)%timeStruct(1)%gru(iGRU)%hru(indxHRU)%var,err,cmessage)
+    end do ! istep
+  end do ! iGRU
+  numGRU = maxGRU-minGRU + 1 
   do iStruct=1,size(structInfo)
     select case(trim(structInfo(iStruct)%structName))
       case('forc')
-        call writeData(ncid,outputTimeStep(indxGRU)%dat(:),outputTimestepUpdate,maxLayers,indxGRU,nSteps,forc_meta, &
-                      outputStructure(1)%forcStat(1),outputStructure(1)%forcStruct(1),'forc', &
+        call writeData(ncid,outputTimeStep(minGRU)%dat(:),outputTimestepUpdate,maxLayers,minGRU,nSteps,&
+                      minGRU, maxGRU, numGRU, & 
+                      forc_meta,outputStructure(1)%forcStat(1),outputStructure(1)%forcStruct(1),'forc', &
                       forcChild_map,outputStructure(1)%indxStruct(1),err,cmessage)
       case('prog')
-        call writeData(ncid,outputTimeStep(indxGRU)%dat(:),outputTimestepUpdate,maxLayers,indxGRU,nSteps,prog_meta, & 
-                      outputStructure(1)%progStat(1),outputStructure(1)%progStruct(1),'prog', &
+        call writeData(ncid,outputTimeStep(minGRU)%dat(:),outputTimestepUpdate,maxLayers,minGRU,nSteps,&
+                      minGRU, maxGRU, numGRU, &
+                      prog_meta,outputStructure(1)%progStat(1),outputStructure(1)%progStruct(1),'prog', &
                       progChild_map,outputStructure(1)%indxStruct(1),err,cmessage)
       case('diag')
-        call writeData(ncid,outputTimeStep(indxGRU)%dat(:),outputTimestepUpdate,maxLayers,indxGRU,nSteps,diag_meta, &
-                      outputStructure(1)%diagStat(1),outputStructure(1)%diagStruct(1),'diag', &
+        call writeData(ncid,outputTimeStep(minGRU)%dat(:),outputTimestepUpdate,maxLayers,minGRU,nSteps,&
+                      minGRU, maxGRU, numGRU, &
+                      diag_meta,outputStructure(1)%diagStat(1),outputStructure(1)%diagStruct(1),'diag', &
                       diagChild_map,outputStructure(1)%indxStruct(1),err,cmessage)
       case('flux')
-        call writeData(ncid,outputTimeStep(indxGRU)%dat(:),outputTimestepUpdate,maxLayers,indxGRU,nSteps,flux_meta, &
-                      outputStructure(1)%fluxStat(1),outputStructure(1)%fluxStruct(1),'flux', &
+        call writeData(ncid,outputTimeStep(minGRU)%dat(:),outputTimestepUpdate,maxLayers,minGRU,nSteps,&
+                      minGRU, maxGRU, numGRU, &
+                      flux_meta,outputStructure(1)%fluxStat(1),outputStructure(1)%fluxStruct(1),'flux', &
                       fluxChild_map,outputStructure(1)%indxStruct(1),err,cmessage)
       case('indx')
-        call writeData(ncid,outputTimeStep(indxGRU)%dat(:),outputTimestepUpdate,maxLayers,indxGRU,nSteps,indx_meta, &
-                      outputStructure(1)%indxStat(1),outputStructure(1)%indxStruct(1),'indx', &
+        call writeData(ncid,outputTimeStep(minGRU)%dat(:),outputTimestepUpdate,maxLayers,minGRU,nSteps,&
+                      minGRU, maxGRU, numGRU, &
+                      indx_meta,outputStructure(1)%indxStat(1),outputStructure(1)%indxStruct(1),'indx', &
                       indxChild_map,outputStructure(1)%indxStruct(1),err,cmessage)
     end select
     if(err/=0)then; message=trim(message)//trim(cmessage)//'['//trim(structInfo(iStruct)%structName)//']'; return; endif
   end do  ! (looping through structures)
   
   do iFreq = 1,maxvarFreq
-    outputTimeStep(indxGRU)%dat(iFreq) = outputTimeStep(indxGRU)%dat(iFreq) + outputTimeStepUpdate(iFreq) 
+    outputTimeStep(minGRU)%dat(iFreq) = outputTimeStep(minGRU)%dat(iFreq) + outputTimeStepUpdate(iFreq) 
   end do ! ifreq
 
 end subroutine
diff --git a/build/source/netcdf/writeOutput.f90 b/build/source/netcdf/writeOutput.f90
index 5733b0f..f1ff324 100644
--- a/build/source/netcdf/writeOutput.f90
+++ b/build/source/netcdf/writeOutput.f90
@@ -150,7 +150,8 @@ end subroutine writeParm
     ! **************************************************************************************
     ! public subroutine writeData: write model time-dependent data
     ! **************************************************************************************
-subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,iGRU,nSteps, & 
+subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,indxGRU,nSteps, &
+            minGRU, maxGRU, numGRU, & 
             meta,stat,dat,structName,map,indx,err,message)
   USE data_types,only:var_info                       ! metadata type
   USE var_lookup,only:maxVarStat                     ! index into stats structure
@@ -159,6 +160,7 @@ subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,iGRU,nSt
   USE var_lookup,only:iLookStat                      ! index into stat structure
   USE globalData,only:outFreq                        ! output file information
   USE globalData,only:outputStructure
+  USE globalData,only:gru_struc
   USE get_ixName_module,only:get_varTypeName         ! to access type strings for error messages
   USE get_ixName_module,only:get_statName            ! to access type strings for error messages
 
@@ -166,10 +168,13 @@ subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,iGRU,nSt
   ! declare dummy variables
   type(var_i)   ,intent(in)        :: ncid              ! file ids
   integer(i4b)  ,intent(inout)     :: outputTimestep(:) ! output time step
-  integer(i4b)  ,intent(inout)        :: outputTimestepUpdate(:) ! number of HRUs in the run domain
+  integer(i4b)  ,intent(inout)     :: outputTimestepUpdate(:) ! number of HRUs in the run domain
   integer(i4b)  ,intent(in)        :: maxLayers         ! maximum number of layers
-  integer(i4b)  ,intent(in)        :: iGRU
-  integer(i4b)  ,intent(in)        :: nSteps             ! number of timeSteps
+  integer(i4b)  ,intent(in)        :: indxGRU
+  integer(i4b)  ,intent(in)        :: nSteps            ! number of timeSteps
+  integer(i4b)  ,intent(in)        :: minGRU            ! minGRU index to write
+  integer(i4b)  ,intent(in)        :: maxGRU            ! maxGRU index to write - probably not needed
+  integer(i4b)  ,intent(in)        :: numGRU            ! number of GRUs to write 
   type(var_info),intent(in)        :: meta(:)           ! meta data
   class(*)      ,intent(in)        :: stat              ! stats data
   class(*)      ,intent(in)        :: dat               ! timestep data
@@ -191,14 +196,16 @@ subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,iGRU,nSt
   integer(i4b)                     :: datLength         ! length of each data vector
   integer(i4b)                     :: maxLength         ! maximum length of each data vector
   real(rkind)                      :: timeVec(nSteps)   ! timeVal to copy
-  real(rkind)                      :: realVec(nSteps)   ! real vector for all HRUs in the run domain
+  real(rkind)                      :: realVec(numGRU, nSteps)   ! real vector for all HRUs in the run domain
   real(rkind)                      :: realArray(nSteps,maxLayers+1)  ! real array for all HRUs in the run domain
   integer(i4b)                     :: intArray(nSteps,maxLayers+1)   ! integer array for all HRUs in the run domain
   integer(i4b)                     :: dataType          ! type of data
   integer(i4b),parameter           :: ixInteger=1001    ! named variable for integer
   integer(i4b),parameter           :: ixReal=1002       ! named variable for real
   integer(i4b)                     :: stepCounter     ! counter to know how much data we have to write
+  integer(i4b)                     :: gruCounter
   integer(i4b)                     :: iStep
+  integer(i4b)                     :: iGRU
   ! initialize error control
   err=0;message="writeData/"
   ! loop through output frequencies
@@ -216,14 +223,10 @@ subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,iGRU,nSt
         call netcdf_err(err,message); if (err/=0) return
         do iStep = 1, nSteps
           ! check if we want this timestep
-          if(.not.outputStructure(1)%finalizeStats(1)%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
+          if(.not.outputStructure(1)%finalizeStats(1)%gru(minGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
           stepCounter = stepCounter+1
-          timeVec(stepCounter) = outputStructure(1)%forcStruct(1)%gru(iGRU)%hru(1)%var(iVar)%tim(iStep)
+          timeVec(stepCounter) = outputStructure(1)%forcStruct(1)%gru(minGRU)%hru(1)%var(iVar)%tim(iStep)
         end do ! iStep
-        ! Write the values
-        ! print*, "Writing Time, startVal", outputTimeStep(iFreq)
-        ! print*, "Writing Time, stepCounter", stepCounter
-        ! print*, "Writing Time, timeVec", timeVec(1:stepCounter)
         err = nf90_put_var(ncid%var(iFreq),ncVarID,timeVec(1:stepCounter),start=(/outputTimestep(iFreq)/),count=(/stepCounter/))
         call netcdf_err(err,message); if (err/=0)then; print*, "err"; return; endif
         ! save the value of the number of steps to update outputTimestep at the end of the function
@@ -236,21 +239,25 @@ subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,iGRU,nSt
       ! check that the variable is desired
       if (iStat==integerMissing.or.trim(meta(iVar)%varName)=='unknown') cycle
 
-        do iHRU=1,gru_struc(iGRU)%hruCount
+        ! do iHRU=1,gru_struc(iGRU)%hruCount
           ! stats output: only scalar variable type
           if(meta(iVar)%varType==iLookVarType%scalarv) then
             select type(stat)
               class is (gru_hru_time_doubleVec)
-                do iStep = 1, nSteps
-                  if(.not.outputStructure(1)%finalizeStats(1)%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
-                  stepCounter = stepCounter + 1
-                  realVec(stepCounter) = stat%gru(iGRU)%hru(iHRU)%var(map(iVar))%tim(iStep)%dat(iFreq)
-                end do
-                ! print*, "MetaData ", meta(iVar)%varName
-                ! print*, "Writing Data, startVal", outputTimeStep(iFreq)
-                ! print*, "Writing Data, stepCounter", stepCounter
-                ! print*, "Writing Data, realVec", realVec(1:stepCounter)
-                err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),realVec(1:stepCounter),start=(/iGRU,outputTimestep(iFreq)/),count=(/1,stepCounter/))
+                
+                gruCounter = 0
+                do iGRU = minGRU, maxGRU
+                  stepCounter = 0
+                  gruCounter = gruCounter + 1
+                  do iStep = 1, nSteps
+                    if(.not.outputStructure(1)%finalizeStats(1)%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
+                    stepCounter = stepCounter + 1
+                    realVec(gruCounter, stepCounter) = stat%gru(iGRU)%hru(1)%var(map(iVar))%tim(iStep)%dat(iFreq)
+                  end do ! iStep
+                  
+                end do ! iGRU
+
+                err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),realVec(1:gruCounter, 1:stepCounter),start=(/minGRU,outputTimestep(iFreq)/),count=(/numGRU,stepCounter/))
                 if (outputTimeStepUpdate(iFreq) /= stepCounter ) then
                   print*, "ERROR Missmatch in Steps"
                   return
@@ -270,9 +277,9 @@ subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,iGRU,nSt
 
 
             ! get the model layers
-            nSoil   = indx%gru(iGRU)%hru(iHRU)%var(iLookIndex%nSoil)%tim(iStep)%dat(1)
-            nSnow   = indx%gru(iGRU)%hru(iHRU)%var(iLookIndex%nSnow)%tim(iStep)%dat(1)
-            nLayers = indx%gru(iGRU)%hru(iHRU)%var(iLookIndex%nLayers)%tim(iStep)%dat(1)
+            nSoil   = indx%gru(iGRU)%hru(1)%var(iLookIndex%nSoil)%tim(iStep)%dat(1)
+            nSnow   = indx%gru(iGRU)%hru(1)%var(iLookIndex%nSnow)%tim(iStep)%dat(1)
+            nLayers = indx%gru(iGRU)%hru(1)%var(iLookIndex%nLayers)%tim(iStep)%dat(1)
 
             ! get the length of each data vector
             select case (meta(iVar)%varType)
@@ -292,14 +299,14 @@ subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,iGRU,nSt
                   do iStep = 1, nSteps
                     if(.not.outputStructure(1)%finalizeStats(1)%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
                     stepCounter = stepCounter + 1
-                    realArray(stepCounter,1:datLength) = dat%gru(iGRU)%hru(iHRU)%var(iVar)%tim(iStep)%dat(:)
+                    realArray(stepCounter,1:datLength) = dat%gru(iGRU)%hru(1)%var(iVar)%tim(iStep)%dat(:)
                   end do
 
                 class is (gru_hru_time_intVec)
                   do iStep = 1, nSteps
                     if(.not.outputStructure(1)%finalizeStats(1)%gru(iGRU)%hru(1)%tim(iStep)%dat(iFreq)) cycle
                     stepCounter = stepCounter + 1
-                    intArray(stepCounter,1:datLength) = dat%gru(iGRU)%hru(iHRU)%var(iVar)%tim(iStep)%dat(:)
+                    intArray(stepCounter,1:datLength) = dat%gru(iGRU)%hru(1)%var(iVar)%tim(iStep)%dat(:)
                   end do
                 class default; err=20; message=trim(message)//'data must not be scalarv and either of type gru_hru_doubleVec or gru_hru_intVec'; return
             end select
@@ -319,20 +326,14 @@ subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,iGRU,nSt
             ! write the data vectors
             select case(dataType)
               case(ixReal)
-                ! print*, "MetaData ", meta(iVar)%varName
-                ! print*, "Writing Data, startVal", outputTimeStep(iFreq)
-                ! print*, "Writing Data, stepCounter", stepCounter
-                ! print*, "Writing Data, realArray", realArray(1:stepCounter,:)
+
                 err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),realArray(1:stepCounter,:),start=(/iGRU,1,outputTimestep(iFreq)/),count=(/1,maxLength,stepCounter/))
                 if (outputTimeStepUpdate(iFreq) /= stepCounter ) then
                   print*, "ERROR Missmatch in Steps"
                   return
                 endif
               case(ixInteger)
-                ! print*, "MetaData ", meta(iVar)%varName
-                ! print*, "Writing Data, startVal", outputTimeStep(iFreq)
-                ! print*, "Writing Data, stepCounter", stepCounter
-                ! print*, "Writing Data, realArray", realArray(1:stepCounter,:)
+
                 err = nf90_put_var(ncid%var(iFreq),meta(iVar)%ncVarID(iFreq),intArray(1:stepCounter,:),start=(/iGRU,1,outputTimestep(iFreq)/),count=(/1,maxLength,stepCounter/))
                 if (outputTimeStepUpdate(iFreq) /= stepCounter ) then
                   print*, "ERROR Missmatch in Steps"
@@ -342,7 +343,7 @@ subroutine writeData(ncid,outputTimestep,outputTimestepUpdate,maxLayers,iGRU,nSt
             end select ! data type
 
           end if ! not scalarv
-        end do ! HRU Loop
+        ! end do ! HRU Loop
 
       ! process error code
       if (err/=0) message=trim(message)//trim(meta(iVar)%varName)//'_'//trim(get_statName(iStat))
-- 
GitLab