diff --git a/build/cmake/build.mac.bash b/build/cmake/build.mac.bash index 78cbb39de..d16d46edc 100755 --- a/build/cmake/build.mac.bash +++ b/build/cmake/build.mac.bash @@ -5,10 +5,13 @@ # Actual settings may vary # Mac Example using MacPorts: -export FC=/opt/local/bin/gfortran # Fortran compiler family +export FC=/opt/homebrew/bin/gfortran # Fortran compiler family #export FLAGS_OPT="-flto=1" # -flto=1 is slow to compile, but might want to use export LIBRARY_LINKS='-llapack' # list of library links -export SUNDIALS_DIR=../../../sundials/instdir/ +export SUNDIALS_DIR=$HOME/local/sundials cmake -B ../cmake_build -S ../. -DUSE_SUNDIALS=ON -DSPECIFY_LAPACK_LINKS=ON -DCMAKE_BUILD_TYPE=Release +#cmake -B ../cmake_build -S ../. -DUSE_SUNDIALS=ON -DSPECIFY_LAPACK_LINKS=ON -DCMAKE_BUILD_TYPE=Debug +#cmake -B ../cmake_build -S ../. -DCMAKE_BUILD_TYPE=Debug \ +# -DCMAKE_Fortran_FLAGS_DEBUG="-O0 -g -fbacktrace -fcheck=all -Wall -Wextra -Wuninitialized -Wmaybe-uninitialized -fcheck=all -finit-real=snan -finit-integer=2147483647" cmake --build ../cmake_build --target all -j diff --git a/build/source/driver/summa_alarms.f90 b/build/source/driver/summa_alarms.f90 index 1ee5d1564..4b3c52d7e 100644 --- a/build/source/driver/summa_alarms.f90 +++ b/build/source/driver/summa_alarms.f90 @@ -21,6 +21,9 @@ module summa_alarms ! used to set alarms to write model output +! named variables for time information +USE globalData, only: numtim ! number of model time steps + ! named variables to define new output files USE globalData, only: noNewFiles ! no new output files USE globalData, only: newFileEveryOct1 ! create a new file on Oct 1 every year (start of the USA water year) @@ -53,19 +56,21 @@ module summa_alarms contains ! used to set alarms to write model output -subroutine summa_setWriteAlarms(oldTime, newTime, endTime, & ! time vectors - newOutputFile, defNewOutputFile, & ! flag to define new output file - ixRestart, printRestart, & ! flag to print the restart file - ixProgress, printProgress, & ! flag to print simulation progress - resetStats, finalizeStats, & ! flags to reset and finalize stats - statCounter, & ! statistics counter - err,message) ! error control +subroutine summa_setWriteAlarms(modelTimeStep, & ! time index + oldTime, newTime, endTime, & ! time vectors + newOutputFile, defNewOutputFile, & ! flag to define new output file + ixRestart, printRestart, & ! flag to print the restart file + ixProgress, printProgress, & ! flag to print simulation progress + resetStats, finalizeStats, & ! flags to reset and finalize stats + statCounter, & ! statistics counter + err,message) ! error control ! --------------------------------------------------------------------------------------- ! data types USE nrtype ! variable types, etc. ! --------------------------------------------------------------------------------------- implicit none ! dummy variables: time vectors + integer(i4b),intent(in) :: modelTimeStep ! index of model time step integer(i4b),intent(in) :: oldTime(:) ! time vector from the previous time step integer(i4b),intent(in) :: newTime(:) ! time vector from the current time step integer(i4b),intent(in) :: endTime(:) ! time vector at the end of the simulation @@ -111,6 +116,12 @@ subroutine summa_setWriteAlarms(oldTime, newTime, endTime, & ! time vector end select + ! check that we do not have multiple files for the buffered write + if(defNewOutputFile .and. modelTimeStep>1)then + err=10 + message=trim(message)//'cannot have multiple output files when using the buffered write decision (check the -n option)'; return + endif + ! ***************************************************************************** ! *** define the need to create a restart file ! ***************************************************************************** @@ -157,6 +168,9 @@ subroutine summa_setWriteAlarms(oldTime, newTime, endTime, & ! time vector case default; err=20; message=trim(message)//'unable to identify output frequency'; return end select + ! force finalize the stats if the last model time step + if(modelTimeStep == numtim) finalizeStats(iFreq)=.true. + ! reset ouput timestep if(resetStats(iFreq)) statCounter(iFreq)=1 @@ -164,4 +178,4 @@ subroutine summa_setWriteAlarms(oldTime, newTime, endTime, & ! time vector end subroutine summa_setWriteAlarms -end module summa_alarms \ No newline at end of file +end module summa_alarms diff --git a/build/source/driver/summa_forcing.f90 b/build/source/driver/summa_forcing.f90 index 79455f302..b94bb0194 100644 --- a/build/source/driver/summa_forcing.f90 +++ b/build/source/driver/summa_forcing.f90 @@ -45,6 +45,8 @@ subroutine summa_readForcing(modelTimeStep, summa1_struc, err, message) ! timing variables USE globalData,only:startRead,endRead ! date/time for the start and end of reading forcing data USE globalData,only:elapsedRead ! elapsed time to read forcing data + ! model decisions + USE globalData,only:model_decisions ! model decision structure ! --------------------------------------------------------------------------------------- ! * variables ! --------------------------------------------------------------------------------------- @@ -73,6 +75,7 @@ subroutine summa_readForcing(modelTimeStep, summa1_struc, err, message) call read_force(& ! input modelTimeStep, & ! intent(in): time step index + model_decisions, & ! intent(in): model decisions structure ! input-output iFile, & ! intent(inout): index of current forcing file in forcing file list forcingStep, & ! intent(inout): index of read position in time dimension in current netcdf file diff --git a/build/source/driver/summa_setup.f90 b/build/source/driver/summa_setup.f90 index 633ca09cf..d8512ae6e 100644 --- a/build/source/driver/summa_setup.f90 +++ b/build/source/driver/summa_setup.f90 @@ -72,7 +72,7 @@ subroutine summa_paramSetup(summa1_struc, err, message) USE nrtype ! variable types, etc. USE summa_type, only:summa1_type_dec ! master summa data type ! subroutines and functions - use time_utils_module,only:elapsedSec ! calculate the elapsed time + USE time_utils_module,only:elapsedSec ! calculate the elapsed time USE mDecisions_module,only:mDecisions ! module to read model decisions USE ffile_info_module,only:ffile_info ! module to read information on forcing datafile USE read_attrb_module,only:read_attrb ! module to read local attributes diff --git a/build/source/driver/summa_writeOutput.f90 b/build/source/driver/summa_writeOutput.f90 index e57309f67..b661cf564 100644 --- a/build/source/driver/summa_writeOutput.f90 +++ b/build/source/driver/summa_writeOutput.f90 @@ -21,10 +21,13 @@ module summa_writeOutput ! used to define/write output files -! named variables to define new output files + ! named variables to define new output files USE globalData, only: noNewFiles ! no new output files USE globalData, only: newFileEveryOct1 ! create a new file on Oct 1 every year (start of the USA water year) +! model decisions +USE globalData,only:model_decisions ! model decision structure + ! metadata USE globalData,only:time_meta ! metadata on the model time USE globalData,only:forc_meta ! metadata on the model forcing data @@ -57,6 +60,47 @@ module summa_writeOutput USE var_lookup,only:iLookPROG ! named variables for local column model prognostic variables USE var_lookup,only:iLookINDEX ! named variables for local column index variables USE var_lookup,only:iLookFREQ ! named variables for the frequency structure +USE var_lookup,only:iLookDECISIONS ! named variables for elements of the decision structure +USE var_lookup,only:iLookVarType ! named variables for variable types + +! generic variable types +USE nrtype ! variable types, etc. + +! data types +USE data_types,only:& + ! final data vectors + dlength, & ! var%dat + ilength, & ! var%dat + ! no spatial dimension + var_i, & ! x%var(:) (i4b) + var_i8, & ! x%var(:) (i8b) + var_d, & ! x%var(:) (dp) + var_ilength, & ! x%var(:)%dat (i4b) + var_dlength, & ! x%var(:)%dat (dp) + ! no variable dimension + hru_i, & ! x%hru(:) (i4b) + hru_d, & ! x%hru(:) (dp) + ! gru dimension + gru_int, & ! x%gru(:)%var(:) (i4b) + gru_double, & ! x%gru(:)%var(:) (dp) + gru_intVec, & ! x%gru(:)%var(:)%dat (i4b) + gru_doubleVec, & ! x%gru(:)%var(:)%dat (dp) + ! gru+hru dimension + gru_hru_int, & ! x%gru(:)%hru(:)%var(:) (i4b) + gru_hru_int8, & ! x%gru(:)%hru(:)%var(:) (i8b) + gru_hru_double, & ! x%gru(:)%hru(:)%var(:) (dp) + gru_hru_intVec, & ! x%gru(:)%hru(:)%var(:)%dat (i4b) + gru_hru_doubleVec ! x%gru(:)%hru(:)%var(:)%dat (dp) + +! metadata structure +USE data_types,only:var_info ! data type for metadata +USE data_types,only:extended_info ! data type for extended metadata +USE data_types,only:struct_info ! summary information on all data structures + +! model write options +USE mDecisions_module,only: & + writePerStep , & ! read forcing data per time step (defualt) + writeFullSeries ! read full forcing series ! safety: set private unless specified otherwise implicit none @@ -70,14 +114,14 @@ subroutine summa_writeOutputFiles(modelTimeStep, summa1_struc, err, message) ! * desired modules ! --------------------------------------------------------------------------------------- ! data types - USE nrtype ! variable types, etc. USE summa_type,only:summa1_type_dec ! master summa data type ! subroutines and functions USE time_utils_module,only:elapsedSec ! calculate the elapsed time + USE mDecisions_module,only:mDecisions ! module to read model decisions USE summa_alarms,only:summa_setWriteAlarms ! set alarms to control model output USE summa_defineOutput,only:summa_defineOutputFiles ! define summa output files USE modelwrite_module,only:writeRestart ! module to write model Restart - USE modelwrite_module,only:writeData,writeBasin ! module to write model output + USE modelwrite_module,only:writeData ! module to write model output USE modelwrite_module,only:writeTime ! module to write model time USE output_stats,only:calcStats ! module for compiling output statistics ! global data: general @@ -91,6 +135,14 @@ subroutine summa_writeOutputFiles(modelTimeStep, summa1_struc, err, message) USE globalData,only:ixProgress ! define frequency to write progress USE globalData,only:ixRestart ! define frequency to write restart files USE globalData,only:newOutputFile ! define option for new output files + ! buffered write + USE globalData,only:numtim ! number of time steps + USE globalData,only:fullIndxSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for indices + USE globalData,only:fullForcSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for forcing + USE globalData,only:fullProgSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for prognostic variables + USE globalData,only:fullDiagSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for diagnostic variables + USE globalData,only:fullFluxSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for flux variables + USE globalData,only:fullBvarSave ! x(:)%gru(:)%var(:) -- saved output for basin variables ! controls on statistics output USE globalData,only:statCounter ! time counter for stats USE globalData,only:resetStats ! flags to reset statistics @@ -118,15 +170,33 @@ subroutine summa_writeOutputFiles(modelTimeStep, summa1_struc, err, message) integer(i4b),intent(out) :: err ! error code character(*),intent(out) :: message ! error message ! local variables - character(LEN=256) :: cmessage ! error message of downwind routine character(len=256) :: timeString ! portion of restart file name that contains the write-out time character(len=256) :: restartFile ! restart file name logical(lgt) :: printRestart=.false. ! flag to print a re-start file logical(lgt) :: printProgress=.false. ! flag to print simulation progress logical(lgt) :: defNewOutputFile=.false. ! flag to define new output files - integer(i4b) :: iGRU,iHRU ! indices of GRUs and HRUs - integer(i4b) :: iStruct ! index of model structure - integer(i4b) :: iFreq ! index of the output frequency + logical(lgt) :: is_writingOutput=.false. ! flag to write model output + logical(lgt) :: is_bufferedWrite=.false. ! flag for buffered write + integer(i4b) :: iTime ! indices of time + integer(i4b) :: iGRU,iHRU ! indices of GRUs and HRUs + integer(i4b) :: iStruct ! index of model structure + integer(i4b) :: iFreq ! index of the output frequency + integer(i4b) :: maxWrite ! maximum number of time steps written + type(gru_hru_int), allocatable :: tempIndxStruct ! Indx temp structure: x%gru(:)%hru(:)%var(:) (rkind) + type(gru_hru_double), allocatable :: tempForcStruct ! Forc temp structure: x%gru(:)%hru(:)%var(:) (rkind) + type(gru_hru_double), allocatable :: tempProgStruct ! Prog temp structure: x%gru(:)%hru(:)%var(:) (rkind) + type(gru_hru_double), allocatable :: tempDiagStruct ! Diag temp structure: x%gru(:)%hru(:)%var(:) (rkind) + type(gru_hru_double), allocatable :: tempFluxStruct ! Flux temp structure: x%gru(:)%hru(:)%var(:) (rkind) + type(gru_double), allocatable :: tempBvarStruct ! Bvar temp structure: x%gru(:)%var(:) (rkind) + type(var_info) , allocatable :: meta(:) ! metadata + type(extended_info) , allocatable :: stat_meta(:) ! statistics metadata (includes only desired variables) + integer(i4b) , allocatable :: child_map(:) ! index of element in child data structure -- meta(map(ivar)) = stat_meta(ivar) + class(*) , allocatable :: timestepData(:) ! vector timestep data (unlimited polymorphic structure) + class(*) , allocatable :: bufferData(:) ! vector buffer data (unlimited polymorphic structure) + class(*) , allocatable :: statsData(:) ! vector stats data (unlimited polymorphic structure) + ! error control + integer(i4b) :: ierr ! error code of downwind routine + character(LEN=256) :: cmessage ! error message of downwind routine ! --------------------------------------------------------------------------------------- ! associate to elements in the data structure summaVars: associate(& @@ -137,7 +207,7 @@ subroutine summa_writeOutputFiles(modelTimeStep, summa1_struc, err, message) diagStat => summa1_struc%diagStat , & ! x%gru(:)%hru(:)%var(:)%dat -- model diagnostic variables fluxStat => summa1_struc%fluxStat , & ! x%gru(:)%hru(:)%var(:)%dat -- model fluxes indxStat => summa1_struc%indxStat , & ! x%gru(:)%hru(:)%var(:)%dat -- model indices - bvarStat => summa1_struc%bvarStat , & ! x%gru(:)%var(:)%dat -- basin-average variabl + bvarStat => summa1_struc%bvarStat , & ! x%gru(:)%var(:)%dat -- basin-average variable ! primary data structures timeStruct => summa1_struc%timeStruct , & ! x%var(:) -- model time data @@ -181,14 +251,43 @@ subroutine summa_writeOutputFiles(modelTimeStep, summa1_struc, err, message) ! initialize number of hru and gru in global data nGRUrun = nGRU nHRUrun = nHRU + endif ! if the first time step + ! ***************************************************************************** + ! *** initialize/populate data structures for the buffered write + ! ***************************************************************************** + + ! if a buffered write + if(model_decisions(iLookDECISIONS%write_buff)%iDecision == writeFullSeries)then + + ! initialize data structures for the buffered write + if(modelTimeStep == 1)then + call initBufferedWrite(structInfo,numtim,ierr,cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); err=20; return; endif + endif ! (if the first time step) + + ! populate data structures + call popBufferStruct(structInfo,summa1_struc,modelTimeStep,ierr,cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); err=20; return; endif + + ! set the maximum number of data points to write + maxWrite = numtim + + else ! (if buffered write) + + ! standard case of write one data value per time step + maxWrite = 1 + + endif ! (if not buffered write) + ! ***************************************************************************** ! *** set alarms for writing data ! ***************************************************************************** ! set alarms to control model output - call summa_setWriteAlarms(oldTime%var, timeStruct%var, finshTime%var, & ! time vectors + call summa_setWriteAlarms(modelTimeStep, & ! time index + oldTime%var, timeStruct%var, finshTime%var, & ! time vectors newOutputFile, defNewOutputFile, & ! flag to define new output file ixRestart, printRestart, & ! flag to print the restart file ixProgress, printProgress, & ! flag to print simulation progress @@ -197,6 +296,17 @@ subroutine summa_writeOutputFiles(modelTimeStep, summa1_struc, err, message) err, cmessage) ! error control if(err/=0)then; message=trim(message)//trim(cmessage); return; endif + ! identify the need to write output file + select case(model_decisions(iLookDECISIONS%write_buff)%iDecision) + case(writePerStep); is_writingOutput = .true. + case(writeFullSeries); is_writingOutput = (modelTimeStep == numtim) + case default + err=10; message=trim(message)//"unknown option for method used to write model output [option="//trim(model_decisions(iLookDECISIONS%write_buff)%cDecision)//"]"; return + end select + + ! check if the buffered write + is_bufferedWrite = (model_decisions(iLookDECISIONS%write_buff)%iDecision == writeFullSeries .and. modelTimeStep == numtim) + ! print progress if(printProgress) write(*,'(i4,1x,5(i2,1x))') timeStruct%var(1:5) @@ -242,36 +352,83 @@ subroutine summa_writeOutputFiles(modelTimeStep, summa1_struc, err, message) call calcStats(bvarStat%gru(iGRU)%var(:),bvarStruct%gru(iGRU)%var(:),statBvar_meta,resetStats,finalizeStats,statCounter,err,cmessage) if(err/=0)then; message=trim(message)//trim(cmessage)//'[bvar stats]'; return; endif - ! write basin-average variables - call writeBasin(iGRU,finalizeStats,outputTimeStep,bvar_meta,bvarStat%gru(iGRU)%var,bvarStruct%gru(iGRU)%var,bvarChild_map,err,cmessage) - if(err/=0)then; message=trim(message)//trim(cmessage)//'[bvar]'; return; endif - end do ! (looping through GRUs) ! **************************************************************************** - ! *** write data + ! *** write integer time information ! **************************************************************************** - ! get the number of HRUs in the run domain - nHRUrun = sum(gru_struc%hruCount) + ! NOTE: This is uncommon -- users rarely require integer time variables (iyyy, im, id, ...) because these can be retrieved from julday - ! write time information - call writeTime(finalizeStats,outputTimeStep,time_meta,timeStruct%var,err,message) + ! NOTE: writing integer time variables is currently restricted to the writePerStep option + if(model_decisions(iLookDECISIONS%write_buff)%iDecision == writePerStep)then + call writeTime(finalizeStats,outputTimeStep,time_meta,timeStruct%var,err,cmessage) + if(err/=0)then; message=trim(message)//trim(cmessage); return; endif + else + if(any(time_meta(:)%varDesire))then + message=trim(message)//"currently no capabilities to output integer time data when using the buffered write option (writeFullSeries)" + err=10; return + endif ! (if desire integer write) + endif ! (if not the writePerStep option) + + ! **************************************************************************** + ! *** write variables for each HRU + ! **************************************************************************** ! write the model output to the NetCDF file - ! Passes the full metadata structure rather than the stats metadata structure because - ! we have the option to write out data of types other than statistics. - ! Thus, we must also pass the stats parent->child maps from childStruct. - do iStruct=1,size(structInfo) - select case(trim(structInfo(iStruct)%structName)) - case('forc'); call writeData(finalizeStats,outputTimeStep,nHRUrun,maxLayers,forc_meta,forcStat,forcStruct,forcChild_map,indxStruct,err,cmessage) - case('prog'); call writeData(finalizeStats,outputTimeStep,nHRUrun,maxLayers,prog_meta,progStat,progStruct,progChild_map,indxStruct,err,cmessage) - case('diag'); call writeData(finalizeStats,outputTimeStep,nHRUrun,maxLayers,diag_meta,diagStat,diagStruct,diagChild_map,indxStruct,err,cmessage) - case('flux'); call writeData(finalizeStats,outputTimeStep,nHRUrun,maxLayers,flux_meta,fluxStat,fluxStruct,fluxChild_map,indxStruct,err,cmessage) - case('indx'); call writeData(finalizeStats,outputTimeStep,nHRUrun,maxLayers,indx_meta,indxStat,indxStruct,indxChild_map,indxStruct,err,cmessage) - end select - if(err/=0)then; message=trim(message)//trim(cmessage)//'['//trim(structInfo(iStruct)%structName)//']'; return; endif - end do ! (looping through structures) + if(is_writingOutput)then + + ! loop through data structures + do iStruct=1,size(structInfo) ! loop means we can apply error code at the end + select case(trim(structInfo(iStruct)%structName)) + + ! define names of desired data structures + case('indx','forc','diag','prog','flux','bvar') ! restrict attention to the variables that we are interested in + + ! get metadata for desired structures + call get_metadata(trim(structInfo(iStruct)%structName), meta, stat_meta, child_map, ierr, cmessage) + if(ierr/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif + + ! get statistics data structure as a 1-element vector (maxWrite=1) + call get_statisticVec(trim(structInfo(iStruct)%structName), maxWrite, summa1_struc, statsData, ierr, cmessage) + if(ierr/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif + + ! ----- write buffered data -------------------------------------------------- + if(is_bufferedWrite)then + + ! get saved data for the buffered write + call get_savedBuffer(trim(structInfo(iStruct)%structName), bufferData, ierr, cmessage) + if(ierr/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif + + ! will only write the (buffered) timestep data but full routine called for convenience + call writeData(is_bufferedWrite,finalizeStats,outputTimeStep,maxWrite,statsData(1),bufferData,meta,child_map,indxStruct,ierr,cmessage) + if(ierr/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif + + ! ----- write data and statistics structures --------------------------------- + else + + ! check per step write + if(maxWrite/=1)then; err=20; message=trim(message)//'expect maxWrite=1'; return; endif + + ! get timestep data structure as a 1-element vector (maxWrite=1) + call get_timestepVec(trim(structInfo(iStruct)%structName), maxWrite, summa1_struc, timestepData, ierr, cmessage) + if(ierr/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif + + ! Passes the full metadata structure and the child map (rather than the stats metadata structure) because + ! we have the option to write out data of types other than statistics. + call writeData(is_bufferedWrite,finalizeStats,outputTimeStep,maxWrite,statsData(1),timestepData,meta,child_map,indxStruct,ierr,cmessage) + if(ierr/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif + + endif ! (write data and statistics structures) + + ! ---------------------------------------------------------------------------- + + ! just keep going if not interested in a data structure + case default; cycle + end select ! select data structure + + end do ! (looping through data structures) + endif ! (if writing output) ! ***************************************************************************** ! *** write restart file @@ -324,6 +481,353 @@ subroutine summa_writeOutputFiles(modelTimeStep, summa1_struc, err, message) end associate summaVars end subroutine summa_writeOutputFiles + + ! ***************************************************************************** + ! ***************************************************************************** + ! ***************************************************************************** + ! ***************************************************************************** + ! ***************************************************************************** + ! ***************************************************************************** + ! private subroutine: initialize data structures for buffered write + subroutine initBufferedWrite(structInfo,numtim,err,message) + ! use modules + USE allocspace_module,only:allocGlobal ! module to allocate space + ! use global data + USE globalData,only:fullIndxSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for indices + USE globalData,only:fullForcSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for forcing + USE globalData,only:fullProgSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for prognostic variables + USE globalData,only:fullDiagSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for diagnostic variables + USE globalData,only:fullFluxSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for flux variables + USE globalData,only:fullBvarSave ! x(:)%gru(:)%var(:) -- saved output for basin variables + implicit none + ! dummy variables + type(struct_info) , intent(in) :: structInfo(:) ! information on the data structures + integer(i4b) , intent(in) :: numtim ! number of model time steps + integer(i4b) , intent(out) :: err ! error code + character(*) , intent(out) :: message ! error message + ! local variables -- temporary data structures + integer(i4b) :: iStruct ! index of data structure + type(gru_hru_int), allocatable :: tempIndx_struct ! Indx temp structure: x%gru(:)%hru(:)%var(:) (i4b) + type(gru_hru_double), allocatable :: tempForc_struct ! Forc temp structure: x%gru(:)%hru(:)%var(:) (rkind) + type(gru_hru_double), allocatable :: tempProg_struct ! Prog temp structure: x%gru(:)%hru(:)%var(:) (rkind) + type(gru_hru_double), allocatable :: tempDiag_struct ! Diag temp structure: x%gru(:)%hru(:)%var(:) (rkind) + type(gru_hru_double), allocatable :: tempFlux_struct ! Flux temp structure: x%gru(:)%hru(:)%var(:) (rkind) + type(gru_double), allocatable :: tempBvar_struct ! Bvar temp structure: x%gru(:)%hru(:)%var(:) (rkind) + ! error control + integer(i4b) :: ierr ! error code of downwind routine + character(LEN=256) :: cmessage ! error message of downwind routine + ! --------------------------------------------------------------------------------------- + ! initialize error control + err=0; message='initBufferedWrite/' + + ! allocate space for local data structures + allocate(tempIndx_struct, tempForc_struct, tempProg_struct, tempDiag_struct, tempFlux_struct, tempBvar_struct, stat=ierr) + if(ierr/=0)then; err=20; message=trim(message)//'problem allocating temporary data structures'; return; endif + + ! allocate space for temporary data structuress + do iStruct=1,size(structInfo) ! loop means we can apply error code at the end + select case(trim(structInfo(iStruct)%structName)) + case('indx'); call allocGlobal(statIndx_meta%var_info, tempIndx_struct, ierr, cmessage) + case('forc'); call allocGlobal(statForc_meta%var_info, tempForc_struct, ierr, cmessage) + case('prog'); call allocGlobal(statProg_meta%var_info, tempProg_struct, ierr, cmessage) + case('diag'); call allocGlobal(statDiag_meta%var_info, tempDiag_struct, ierr, cmessage) + case('flux'); call allocGlobal(statFlux_meta%var_info, tempFlux_struct, ierr, cmessage) + case('bvar'); call allocGlobal(statBvar_meta%var_info, tempBvar_struct, ierr, cmessage) + case default; cycle ! do not expect additional data structures + end select + if(ierr/=0)then; message=trim(message)//trim(cmessage)//'['//trim(structInfo(iStruct)%structName)//']'; err=20; return; endif + end do ! (looping through structures) + + ! add a time dimension + do iStruct=1,size(structInfo) ! loop means we can apply error code at the end + select case(trim(structInfo(iStruct)%structName)) + case('indx'); allocate(fullIndxSave(numtim), source=tempIndx_struct, stat=ierr) + case('forc'); allocate(fullForcSave(numtim), source=tempForc_struct, stat=ierr) + case('prog'); allocate(fullProgSave(numtim), source=tempProg_struct, stat=ierr) + case('diag'); allocate(fullDiagSave(numtim), source=tempDiag_struct, stat=ierr) + case('flux'); allocate(fullFluxSave(numtim), source=tempFlux_struct, stat=ierr) + case('bvar'); allocate(fullBvarSave(numtim), source=tempBvar_struct, stat=ierr) + case default; cycle ! do not expect additional data structures + end select + if(ierr/=0)then; message=trim(message)//trim(cmessage)//'['//trim(structInfo(iStruct)%structName)//']'; err=20; return; endif + end do ! (looping through structues) + + ! deallocate space for data structures + deallocate(tempIndx_struct, tempForc_struct, tempProg_struct, tempDiag_struct, tempFlux_struct, tempBvar_struct, stat=ierr) + if(ierr/=0)then; err=20; message=trim(message)//'problem deallocating temporary data structures'; return; endif + + end subroutine initBufferedWrite + + ! ***************************************************************************** + ! ***************************************************************************** + ! private subroutine: populate data structures for buffered write + subroutine popBufferStruct(structInfo,summaStruct,iTime,err,message) + ! global data: structures for buffered write + USE globalData,only:fullIndxSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for indices + USE globalData,only:fullForcSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for forcing + USE globalData,only:fullProgSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for prognostic variables + USE globalData,only:fullDiagSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for diagnostic variables + USE globalData,only:fullFluxSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for flux variables + USE globalData,only:fullBvarSave ! x(:)%gru(:)%var(:) -- saved output for basin variables + ! global data: structures for GRU-HRU topology + USE globalData,only:gru_struc ! gru-hru mapping structures + ! derived types + USE summa_type,only:summa1_type_dec ! master summa data type + implicit none + ! input variables + type(struct_info) , intent(in) :: structInfo(:) ! information on the data structures + type(summa1_type_dec), intent(in) :: summaStruct ! master summa data structure + integer(i4b) , intent(in) :: iTime ! index of time (modelTimeStep) + ! output variables + integer(i4b) , intent(out) :: err ! error code + character(*) , intent(out) :: message ! error message + ! local variables + integer(i4b) :: iStruct ! index of data structure + integer(i4b) :: iGRU ! index of GRU + integer(i4b) :: iHRU ! index of HRU + integer(i4b) :: iVar ! index of variable + integer(i4b) :: pVar ! index of "parent" variable (i.e., index in the data structure) + type(var_info) , allocatable :: meta(:) ! metadata + type(extended_info) , allocatable :: stat_meta(:) ! statistics metadata (includes only desired variables) + integer(i4b) , allocatable :: child_map(:) ! index of element in child data structure -- meta(map(ivar)) = stat_meta(ivar) + ! error control + integer(i4b) :: ierr ! local error code + character(len=256) :: cmessage ! error message of the downwind routine + ! associate to elements in the data structure + ! ---------------------------------------------------------------------------------------------------------------------------- + ! primary data structures + summaAssociate: associate(& + indxStruct => summaStruct%indxStruct , & ! x%gru(:)%hru(:)%var(:)%dat -- model indices + forcStruct => summaStruct%forcStruct , & ! x%gru(:)%hru(:)%var(:) -- model forcing data + progStruct => summaStruct%progStruct , & ! x%gru(:)%hru(:)%var(:)%dat -- model prognostic (state) variables + diagStruct => summaStruct%diagStruct , & ! x%gru(:)%hru(:)%var(:)%dat -- model diagnostic variables + fluxStruct => summaStruct%fluxStruct , & ! x%gru(:)%hru(:)%var(:)%dat -- model fluxes + bvarStruct => summaStruct%bvarStruct , & ! x%gru(:)%var(:)%dat -- basin-average variables + nGRU => summaStruct%nGRU & + ) ! assignment to variables in the data structures + ! ------------------------------------------------------------------------------------------------------------------------- + ! initialize error control + err=0; message='popBufferStruct/' + + ! loop through GRUs and HRUs + do iGRU=1,nGRU + do iHRU=1,gru_struc(iGRU)%hruCount + + ! loop through data structures + do iStruct=1,size(structInfo) ! loop means we can apply error code at the end + + ! ----- get metadata for desired structure ----------------------------------------------------------------------------- + + ! select structure + select case(trim(structInfo(iStruct)%structName)) + + ! get metadata for desired structures + case('indx','forc','diag','prog','flux','bvar') ! restrict attention to the variables that we are interested in + call get_metadata(trim(structInfo(iStruct)%structName), meta, stat_meta, child_map, ierr, cmessage) + if(ierr>0)then; err=20; message=trim(message)//trim(cmessage); return; endif + + ! just keep going if not interested in a data structure + case default; cycle + + end select ! select data structure + + ! ----- populate data for desired variables ---------------------------------------------------------------------------- + + ! NOTE: use stat_meta because these are the variables desired in the output + do iVar=1,size(stat_meta) ! skip if size=0 + + ! don't do anything if var is not requested or not a scalar variable + if (.not.stat_meta(iVar)%varDesire) cycle ! variable not requested + if (stat_meta(iVar)%varType/=iLookVarType%outstat) cycle ! not a scalar variable + + ! index in parent structure + pVar = stat_meta(iVar)%ixParent + if(trim(stat_meta(iVar)%varName) /= trim(meta(pVar)%varName))then + message=trim(message)//'variable names do not match' + err=20; return + endif + + ! populate GRU+HRU structures + select case(trim(structInfo(iStruct)%structName)) + case('indx'); fullIndxSave(iTime)%gru(iGRU)%hru(iHRU)%var(iVar) = indxStruct%gru(iGRU)%hru(iHRU)%var(pVar)%dat(1) + case('forc'); fullForcSave(iTime)%gru(iGRU)%hru(iHRU)%var(iVar) = forcStruct%gru(iGRU)%hru(iHRU)%var(pVar) + case('diag'); fullProgSave(iTime)%gru(iGRU)%hru(iHRU)%var(iVar) = progStruct%gru(iGRU)%hru(iHRU)%var(pVar)%dat(1) + case('prog'); fullDiagSave(iTime)%gru(iGRU)%hru(iHRU)%var(iVar) = diagStruct%gru(iGRU)%hru(iHRU)%var(pVar)%dat(1) + case('flux'); fullFluxSave(iTime)%gru(iGRU)%hru(iHRU)%var(iVar) = fluxStruct%gru(iGRU)%hru(iHRU)%var(pVar)%dat(1) + case('bvar') ! GRU-only data structure + if(iHRU==1) fullBvarSave(iTime)%gru(iGRU)%var(iVar) = bvarStruct%gru(iGRU)%var(pVar)%dat(1) + case default; err=20; message=trim(message)//'do not expect any other structures than what is listed'; return + end select + + end do ! (looping through variables) + + end do ! (looping through structures) + + ! deallocate space for metadata structures + deallocate(meta, stat_meta, child_map, stat=ierr) + if(ierr/=0)then; message=trim(message)//'problem allocating metadata'; err=20; return; endif + + end do ! (looping through HRUs) + end do ! (looping through GRUs) + + end associate summaAssociate + + end subroutine popBufferStruct + + ! ***************************************************************************** + ! ***************************************************************************** + + subroutine get_metadata(tag, meta, stat_meta, child_map, err, message) + ! dummy variables + character(*) , intent(in) :: tag ! name of data structure + type(var_info) , intent(out), allocatable :: meta(:) ! metadata + type(extended_info) , intent(out), allocatable :: stat_meta(:) ! statistics metadata + integer(i4b) , intent(out), allocatable :: child_map(:) ! index of the child data structure + integer(i4b) , intent(out) :: err ! error code + character(*) , intent(out) :: message ! error message + ! local variables + integer(i4b) :: nVar ! number of variables in the metadata structure + integer(i4b) :: nSubset ! number of variables in the statistics metadata structure + integer(i4b) :: ierr ! local error code + ! initalize error control + err=0; message='get_metadata/' + + ! get size of metadata structures for a given structure + select case(trim(tag)) + case('indx'); nVar = size(indx_meta); nSubset = size(statIndx_meta) + case('forc'); nVar = size(forc_meta); nSubset = size(statForc_meta) + case('prog'); nVar = size(prog_meta); nSubset = size(statProg_meta) + case('diag'); nVar = size(diag_meta); nSubset = size(statDiag_meta) + case('flux'); nVar = size(flux_meta); nSubset = size(statFlux_meta) + case('bvar'); nVar = size(bvar_meta); nSubset = size(statBvar_meta) + case default; err=10; message=trim(message)//'data structure not needed -- should not get here'; return + end select + + ! allocate space for the metadata structures + allocate(meta(nVar), stat_meta(nSubset), child_map(nVar), stat=ierr) + if(ierr/=0)then; message=trim(message)//'problem allocating metadata'; err=20; return; endif + + ! save metadata + select case(trim(tag)) + case('indx'); meta(:) = indx_meta(:); stat_meta(:) = statIndx_meta(:); child_map(:) = indxChild_map(:) + case('forc'); meta(:) = forc_meta(:); stat_meta(:) = statForc_meta(:); child_map(:) = forcChild_map(:) + case('prog'); meta(:) = prog_meta(:); stat_meta(:) = statProg_meta(:); child_map(:) = progChild_map(:) + case('diag'); meta(:) = diag_meta(:); stat_meta(:) = statDiag_meta(:); child_map(:) = diagChild_map(:) + case('flux'); meta(:) = flux_meta(:); stat_meta(:) = statFlux_meta(:); child_map(:) = fluxChild_map(:) + case('bvar'); meta(:) = bvar_meta(:); stat_meta(:) = statBvar_meta(:); child_map(:) = bvarChild_map(:) + end select + + end subroutine get_metadata + + ! ***************************************************************************** + ! ***************************************************************************** + + subroutine get_savedBuffer(tag, structVec, err, message) + ! saved timestep data + USE globalData,only:fullIndxSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for indices + USE globalData,only:fullForcSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for forcing + USE globalData,only:fullProgSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for prognostic variables + USE globalData,only:fullDiagSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for diagnostic variables + USE globalData,only:fullFluxSave ! x(:)%gru(:)%hru(:)%var(:) -- saved output for flux variables + USE globalData,only:fullBvarSave ! x(:)%gru(:)%var(:) -- saved output for basin variables + ! dummy variables + character(*) , intent(in) :: tag ! name of data structure + class(*) , intent(out), allocatable :: structVec(:) ! vector of any data structure + integer(i4b) , intent(out) :: err ! error code + character(*) , intent(out) :: message ! error message + ! initialize errror control + err=0; message='get_savedBuffer/' + + ! allocate data + select case(trim(tag)) + case('indx'); allocate(structVec, source=fullIndxSave, stat=err) + case('forc'); allocate(structVec, source=fullForcSave, stat=err) + case('prog'); allocate(structVec, source=fullProgSave, stat=err) + case('diag'); allocate(structVec, source=fullDiagSave, stat=err) + case('flux'); allocate(structVec, source=fullFluxSave, stat=err) + case('bvar'); allocate(structVec, source=fullBvarSave, stat=err) + case default; err=20; message=trim(message)//'cannot identify data structure'; return + end select + + ! check errors + if(err/=0)then + message=trim(message)//'problem allocating space for structure '//trim(tag) + err=20; return + endif + + end subroutine get_savedBuffer + + ! ***************************************************************************** + ! ***************************************************************************** + + subroutine get_timestepVec(tag, nVec, summaStruct, structVec, err, message) + USE summa_type,only:summa1_type_dec ! master summa data type + ! dummy variables + character(*) , intent(in) :: tag ! name of data structure + integer(i4b) , intent(in) :: nVec ! number of vector elements + type(summa1_type_dec) , intent(in) :: summaStruct ! master summa data structure + class(*) , intent(out), allocatable :: structVec(:) ! vector of any data structure + integer(i4b) , intent(out) :: err ! error code + character(*) , intent(out) :: message ! error message + ! initialize errror control + err=0; message='get_timestepVec/' + + ! allocate data + select case(trim(tag)) + case('indx'); allocate(structVec(nVec), source=summaStruct%indxStruct, stat=err) + case('forc'); allocate(structVec(nVec), source=summaStruct%forcStruct, stat=err) + case('prog'); allocate(structVec(nVec), source=summaStruct%progStruct, stat=err) + case('diag'); allocate(structVec(nVec), source=summaStruct%diagStruct, stat=err) + case('flux'); allocate(structVec(nVec), source=summaStruct%fluxStruct, stat=err) + case('bvar'); allocate(structVec(nVec), source=summaStruct%bvarStruct, stat=err) + case default; err=20; message=trim(message)//'cannot identify data structure'; return + end select + + ! check errors + if(err/=0)then + message=trim(message)//'problem allocating space for structure '//trim(tag) + err=20; return + endif + + end subroutine get_timestepVec + + ! ***************************************************************************** + ! ***************************************************************************** + + subroutine get_statisticVec(tag, nVec, summaStruct, structVec, err, message) + USE summa_type,only:summa1_type_dec ! master summa data type + ! dummy variables + character(*) , intent(in) :: tag ! name of data structure + integer(i4b) , intent(in) :: nVec ! number of vector elements + type(summa1_type_dec) , intent(in) :: summaStruct ! master summa data structure + class(*) , intent(out), allocatable :: structVec(:) ! vector of any data structure + integer(i4b) , intent(out) :: err ! error code + character(*) , intent(out) :: message ! error message + ! initialize errror control + err=0; message='get_statisticVec/' + + ! allocate data + select case(trim(tag)) + case('indx'); allocate(structVec(nVec), source=summaStruct%indxStat, stat=err) + case('forc'); allocate(structVec(nVec), source=summaStruct%forcStat, stat=err) + case('prog'); allocate(structVec(nVec), source=summaStruct%progStat, stat=err) + case('diag'); allocate(structVec(nVec), source=summaStruct%diagStat, stat=err) + case('flux'); allocate(structVec(nVec), source=summaStruct%fluxStat, stat=err) + case('bvar'); allocate(structVec(nVec), source=summaStruct%bvarStat, stat=err) + case default; err=20; message=trim(message)//'cannot identify data structure'; return + end select + + ! check errors + if(err/=0)then + message=trim(message)//'problem allocating space for structure '//trim(tag) + err=20; return + endif + + end subroutine get_statisticVec + + ! ***************************************************************************** + ! ***************************************************************************** + end module summa_writeOutput diff --git a/build/source/dshare/get_ixname.f90 b/build/source/dshare/get_ixname.f90 index 4c3fddd1a..71b602b53 100644 --- a/build/source/dshare/get_ixname.f90 +++ b/build/source/dshare/get_ixname.f90 @@ -100,6 +100,8 @@ function get_ixdecisions(varName) case('infRateMax' ); get_ixdecisions=iLookDECISIONS%infRateMax ! choice of maximum infiltration rate method case('surfRun_IE' ); get_ixdecisions=iLookDECISIONS%surfRun_IE ! choice of parameterization for infiltration excess surface runoff case('surfRun_SE' ); get_ixdecisions=iLookDECISIONS%surfRun_SE ! choice of parameterization for saturation excess surface runoff + case('read_force' ); get_ixdecisions=iLookDECISIONS%read_force ! method used to read forcing data (per step or full read) + case('write_buff' ); get_ixdecisions=iLookDECISIONS%write_buff ! method used to buffer writing of model output (none, full) ! get to here if cannot find the variable case default get_ixdecisions = integerMissing @@ -1130,7 +1132,7 @@ subroutine get_ixUnknown(varName,typeName,vDex,err,message) end subroutine get_ixUnknown ! ******************************************************************************************************************* - ! public function get_ixFreq: get the index of the named variables for the output frequencies + ! public function get_ixLookup: get the index of the named variables for lookup ! ******************************************************************************************************************* function get_ixLookup(varName) USE var_lookup,only:iLookLOOKUP ! indices of the named variables diff --git a/build/source/dshare/globalData.f90 b/build/source/dshare/globalData.f90 index 812c428ad..a8a368e7e 100644 --- a/build/source/dshare/globalData.f90 +++ b/build/source/dshare/globalData.f90 @@ -37,6 +37,9 @@ MODULE globalData USE data_types,only:extended_info ! extended metadata for variables in each model structure USE data_types,only:struct_info ! summary information on all data structures USE data_types,only:var_i ! vector of integers + USE data_types,only:gru_hru_int ! x%gru(:)%hru(:)%var(:) (i4b) + USE data_types,only:gru_hru_double ! x%gru(:)%hru(:)%var(:) (rkind) + USE data_types,only:gru_double ! x%gru(:)%var(:) (rkind) ! number of variables in each data structure USE var_lookup,only:maxvarTime ! time: maximum number variables USE var_lookup,only:maxvarForc ! forcing data: maximum number variables @@ -232,6 +235,17 @@ MODULE globalData integer(i4b),save,public :: chunksize=1024 ! chunk size for the netcdf read/write integer(i4b),save,public :: outputPrecision=nf90_double ! variable type integer(i4b),save,public :: outputCompressionLevel=4 ! output netcdf file deflate level: 0-9. 0 is no compression. + ! define data structures for the buffered read + integer(i4b),save,public :: ixStartRead ! start index of the data read + real(rkind),save,public,allocatable :: fulltimeVec(:) ! full time vector in an input file (nRead) + type(gru_hru_double),save,public,allocatable :: fullforcingStruct(:) ! x(:)%gru(:)%hru(:)%var(:) -- full model forcing data + ! define data structures for the buffered write + type(gru_hru_int), save,public,allocatable :: fullIndxSave(:) ! x(:)%gru(:)%hru(:)%var(:) -- saved output for indices + type(gru_hru_double),save,public,allocatable :: fullForcSave(:) ! x(:)%gru(:)%hru(:)%var(:) -- saved output for forcing + type(gru_hru_double),save,public,allocatable :: fullProgSave(:) ! x(:)%gru(:)%hru(:)%var(:) -- saved output for prognostic variables + type(gru_hru_double),save,public,allocatable :: fullDiagSave(:) ! x(:)%gru(:)%hru(:)%var(:) -- saved output for diagnostic variables + type(gru_hru_double),save,public,allocatable :: fullFluxSave(:) ! x(:)%gru(:)%hru(:)%var(:) -- saved output for flux variables + type(gru_double), save,public,allocatable :: fullBvarSave(:) ! x(:)%gru(:)%var(:) -- saved output for basin variables ! define result from the time calls integer(i4b),dimension(8),save,public :: startInit,endInit ! date/time for the start and end of the initialization integer(i4b),dimension(8),save,public :: startSetup,endSetup ! date/time for the start and end of the parameter setup @@ -239,7 +253,7 @@ MODULE globalData integer(i4b),dimension(8),save,public :: startRead,endRead ! date/time for the start and end of the data read integer(i4b),dimension(8),save,public :: startWrite,endWrite ! date/time for the start and end of the stats/write integer(i4b),dimension(8),save,public :: startPhysics,endPhysics ! date/time for the start and end of the physics - ! define elapsed time + ! define elapsed time real(rkind),save,public :: elapsedInit ! elapsed time for the initialization real(rkind),save,public :: elapsedSetup ! elapsed time for the parameter setup real(rkind),save,public :: elapsedRestart ! elapsed time to read restart data diff --git a/build/source/dshare/outpt_stat.f90 b/build/source/dshare/outpt_stat.f90 index 56fafbd3a..e67c165e8 100644 --- a/build/source/dshare/outpt_stat.f90 +++ b/build/source/dshare/outpt_stat.f90 @@ -54,7 +54,7 @@ subroutine calcStats(stat,dat,meta,resetStats,finalizeStats,statCounter,err,mess character(256) :: cmessage ! error message integer(i4b) :: iVar ! index for varaiable loop integer(i4b) :: pVar ! index into parent structure - real(rkind) :: tdata ! dummy for pulling info from dat structure + real(rkind) :: tdata ! dummy for pulling info from dat structure ! initialize error control err=0; message='calcStats/' @@ -73,9 +73,9 @@ subroutine calcStats(stat,dat,meta,resetStats,finalizeStats,statCounter,err,mess ! extract data from the structures select type (dat) - type is (real(rkind)); tdata = dat(pVar) - class is (dlength) ; tdata = dat(pVar)%dat(1) - class is (ilength) ; tdata = real(dat(pVar)%dat(1), kind(rkind)) + type is (real(rkind)); tdata = dat(pVar) + class is (dlength) ; tdata = dat(pVar)%dat(1) + class is (ilength) ; tdata = real(dat(pVar)%dat(1), kind(rkind)) class default;err=20;message=trim(message)//'dat type not found';return end select @@ -114,7 +114,7 @@ subroutine calc_stats(meta,stat,tdata,resetStats,finalizeStats,statCounter,err,m ! input variables class(var_info),intent(in) :: meta ! meta data structure class(*) ,intent(inout) :: stat ! statistics structure - real(rkind) ,intent(in) :: tdata ! data value + real(rkind) ,intent(in) :: tdata ! data value logical(lgt) ,intent(in) :: resetStats(:) ! vector of flags to reset statistics logical(lgt) ,intent(in) :: finalizeStats(:) ! vector of flags to reset statistics integer(i4b) ,intent(in) :: statCounter(:) ! number of time steps in each output frequency @@ -122,8 +122,8 @@ subroutine calc_stats(meta,stat,tdata,resetStats,finalizeStats,statCounter,err,m integer(i4b) ,intent(out) :: err ! error code character(*) ,intent(out) :: message ! error message ! internals - real(rkind),dimension(maxvarFreq*2) :: tstat ! temporary stats vector - integer(i4b) :: iFreq ! index of output frequency + real(rkind),dimension(maxvarFreq*2) :: tstat ! temporary stats vector + integer(i4b) :: iFreq ! index of output frequency ! initialize error control err=0; message='calc_stats/' diff --git a/build/source/dshare/popMetadat.f90 b/build/source/dshare/popMetadat.f90 index 391833ae5..85417b149 100644 --- a/build/source/dshare/popMetadat.f90 +++ b/build/source/dshare/popMetadat.f90 @@ -1008,6 +1008,7 @@ subroutine read_output_file(err,message) ! option 0: file format = varName | outFreq ! option 1: file format = varName | outFreq | statisticName ! option 2: file format = varName | outFreq | inst | sum | mean | var | min | max | mode + ! e.g., varName | outFreq | 0 | 0 | 1 | 0 | 0 | 0 | 0 select case(nWords) case(nameIndex + 2, nameIndex); fileFormat=noStatsDesired ! no statistic desired (temporally constant variables) case(freqIndex + 2 ); fileFormat=provideStatName ! provide the name of the desired statistic diff --git a/build/source/dshare/var_lookup.f90 b/build/source/dshare/var_lookup.f90 index d2ec2d03e..1c495d061 100644 --- a/build/source/dshare/var_lookup.f90 +++ b/build/source/dshare/var_lookup.f90 @@ -79,6 +79,8 @@ MODULE var_lookup integer(i4b) :: infRateMax = integerMissing ! choice of method to determine maximum infiltration rate integer(i4b) :: surfRun_IE = integerMissing ! choice of parameterization for infiltration excess surface runoff integer(i4b) :: surfRun_SE = integerMissing ! choice of parameterization for saturation excess surface runoff + integer(i4b) :: read_force = integerMissing ! method used to read forcing data (per step or full read) + integer(i4b) :: write_buff = integerMissing ! method used to buffer model write (none, per file) endtype iLook_decision @@ -906,7 +908,7 @@ MODULE var_lookup 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,& 21, 22, 23, 24, 25, 26, 27, 28, 29, 30,& 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,& - 41, 42, 43) + 41, 42, 43, 44, 45) ! named variables: model time type(iLook_time), public,parameter :: iLookTIME =iLook_time ( 1, 2, 3, 4, 5, 6, 7) ! named variables: model forcing data diff --git a/build/source/engine/allocspace.f90 b/build/source/engine/allocspace.f90 index dd77f4eb9..0f421903e 100644 --- a/build/source/engine/allocspace.f90 +++ b/build/source/engine/allocspace.f90 @@ -117,7 +117,7 @@ subroutine allocGlobal(metaStruct,dataStruct,err,message) class is (gru_hru_intVec); if(allocated(dataStruct%gru))then; check=.true.; else; allocate(dataStruct%gru(nGRU),stat=err); end if class is (gru_hru_double); if(allocated(dataStruct%gru))then; check=.true.; else; allocate(dataStruct%gru(nGRU),stat=err); end if class is (gru_hru_doubleVec); if(allocated(dataStruct%gru))then; check=.true.; else; allocate(dataStruct%gru(nGRU),stat=err); end if - ! gru+hru+z dimensions + ! gru+hru+z dimensions class is (gru_hru_z_vLookup); if(allocated(dataStruct%gru))then; check=.true.; else; allocate(dataStruct%gru(nGRU),stat=err); end if end select @@ -134,7 +134,7 @@ subroutine allocGlobal(metaStruct,dataStruct,err,message) class is (gru_hru_intVec); if(allocated(dataStruct%gru(iGRU)%hru))then; check=.true.; else; allocate(dataStruct%gru(iGRU)%hru(gru_struc(iGRU)%hruCount),stat=err); end if class is (gru_hru_double); if(allocated(dataStruct%gru(iGRU)%hru))then; check=.true.; else; allocate(dataStruct%gru(iGRU)%hru(gru_struc(iGRU)%hruCount),stat=err); end if class is (gru_hru_doubleVec); if(allocated(dataStruct%gru(iGRU)%hru))then; check=.true.; else; allocate(dataStruct%gru(iGRU)%hru(gru_struc(iGRU)%hruCount),stat=err); end if - class is (gru_hru_z_vLookup); if(allocated(dataStruct%gru(iGRU)%hru))then; check=.true.; else; allocate(dataStruct%gru(iGRU)%hru(gru_struc(iGRU)%hruCount),stat=err); end if + class is (gru_hru_z_vLookup); if(allocated(dataStruct%gru(iGRU)%hru))then; check=.true.; else; allocate(dataStruct%gru(iGRU)%hru(gru_struc(iGRU)%hruCount),stat=err); end if class default ! do nothing: It is acceptable to not be any of these specified cases end select ! check errors @@ -533,7 +533,7 @@ end subroutine copyStruct_i4b ! private subroutine allocateDat_rkind: initialize data dimension of the data structures ! ************************************************************************************************ subroutine allocateDat_rkind(metadata,nSnow,nSoil,nLayers, & ! input - varData,err,message) ! output + varData,err,message) ! output ! access subroutines USE get_ixName_module,only:get_varTypeName ! to access type strings for error messages diff --git a/build/source/engine/mDecisions.f90 b/build/source/engine/mDecisions.f90 index 517a8f051..52bb5d1d6 100644 --- a/build/source/engine/mDecisions.f90 +++ b/build/source/engine/mDecisions.f90 @@ -167,6 +167,12 @@ module mDecisions_module integer(i4b),parameter,public :: FUSEPRMS = 353 ! FUSE PRMS surface runoff integer(i4b),parameter,public :: FUSEAVIC = 354 ! FUSE ARNO/VIC surface runoff integer(i4b),parameter,public :: FUSETOPM = 355 ! FUSE TOPMODEL surface runoff +! look-up values for the buffered read of forcing data +integer(i4b),parameter,public :: readPerStep = 361 ! read forcing data per time step (defualt) +integer(i4b),parameter,public :: readFullSeries = 362 ! read full forcing series +! look-up values for the buffered write of model output +integer(i4b),parameter,public :: writePerStep = 371 ! write data per time step (defualt) +integer(i4b),parameter,public :: writeFullSeries = 372 ! write all data for a given output file ! ----------------------------------------------------------------------------------------------------------- @@ -713,6 +719,24 @@ subroutine mDecisions(err,message) err=10; message=trim(message)//"unknown option for saturation excess surface runoff method [option="//trim(model_decisions(iLookDECISIONS%surfRun_SE)%cDecision)//"]"; return end select + ! method used to read forcing data (per step or full read) + ! NOTE: use read forcing data per time step as the defualt + select case(trim(model_decisions(iLookDECISIONS%read_force)%cDecision)) + case('readPerStep','notPopulatedYet'); model_decisions(iLookDECISIONS%read_force)%iDecision = readPerStep ! read forcing data per time step (defualt) + case('readFullSeries' ); model_decisions(iLookDECISIONS%read_force)%iDecision = readFullSeries ! read full forcing series + case default + err=10; message=trim(message)//"unknown option for method used to read forcing data [option="//trim(model_decisions(iLookDECISIONS%read_force)%cDecision)//"]"; return + end select + + ! method used to write model output (per step or full write) + ! NOTE: use per time step as the defualt + select case(trim(model_decisions(iLookDECISIONS%write_buff)%cDecision)) + case('writePerStep','notPopulatedYet'); model_decisions(iLookDECISIONS%write_buff)%iDecision = writePerStep ! write model output per time step (defualt) + case('writeFullSeries' ); model_decisions(iLookDECISIONS%write_buff)%iDecision = writeFullSeries ! write all data for a given output file + case default + err=10; message=trim(message)//"unknown option for method used to write model output [option="//trim(model_decisions(iLookDECISIONS%write_buff)%cDecision)//"]"; return + end select + ! ----------------------------------------------------------------------------------------------------------------------------------------------- ! check for consistency among options ! ----------------------------------------------------------------------------------------------------------------------------------------------- diff --git a/build/source/engine/read_force.f90 b/build/source/engine/read_force.f90 index f6a34335d..a4940930b 100644 --- a/build/source/engine/read_force.f90 +++ b/build/source/engine/read_force.f90 @@ -24,7 +24,9 @@ module read_force_module USE nrtype ! variable types, etc. ! derived data types +USE data_types,only:var_ilength ! x%var(:)%dat(:) (i4b) USE data_types,only:gru_hru_double ! x%gru(:)%hru(:)%var(:) (rkind) +USE data_types,only:model_options ! defines the model decisions ! constants USE multiconst,only:secprday ! number of seconds in a day @@ -40,6 +42,7 @@ module read_force_module USE globalData,only:ixHRUfile_min,ixHRUfile_max ! global data on the forcing file +USE globalData,only:numtim ! number time steps USE globalData,only:data_step ! length of the data step (s) USE globalData,only:forcFileInfo ! forcing file info USE globalData,only:dJulianStart ! julian day of start time of simulation @@ -49,11 +52,21 @@ module read_force_module USE globalData,only:yearLength ! number of days in the current year USE globalData,only:nHRUfile ! number of days in the data file +! global data holding the buffered read +USE globalData,only:ixStartRead ! start index of the data read +USE globalData,only:fulltimeVec ! full time vector in an input file (nRead) +USE globalData,only:fullforcingStruct ! full forcing data structure + ! global metadata USE globalData,only:time_meta,forc_meta ! metadata structures USE var_lookup,only:iLookTIME,iLookFORCE ! named variables to define structure elements USE var_lookup,only:iLookDECISIONS ! named variables for elements of the decision structure +! look-up values for option to read input data +USE mDecisions_module,only: & + readPerStep , & ! read forcing data per time step (defualt) + readFullSeries ! read full forcing series + ! file paths USE summaFileManager,only:FORCING_PATH ! path of the forcing data file @@ -68,39 +81,54 @@ module read_force_module contains - ! ************************************************************************************************ ! public subroutine read_force: read in forcing data ! ************************************************************************************************ - subroutine read_force(istep,iFile,iRead,ncid,time_data,forcStruct,err,message) + subroutine read_force(istep,model_decisions,iFile,iRead,ncid,time_data,forcStruct,err,message) ! provide access to subroutines - USE netcdf ! netcdf capability - USE time_utils_module,only:compJulDay ! convert calendar date to julian day - USE time_utils_module,only:compcalday ! convert julian day to calendar date - USE time_utils_module,only:elapsedSec ! calculate the elapsed time + USE netcdf ! netcdf capability + USE time_utils_module,only:compJulDay ! convert calendar date to julian day + USE time_utils_module,only:compcalday ! convert julian day to calendar date + USE time_utils_module,only:elapsedSec ! calculate the elapsed time implicit none ! define input variables - integer(i4b),intent(in) :: istep ! time index AFTER the start index + integer(i4b),intent(in) :: istep ! time index AFTER the start index + type(model_options),intent(in) :: model_decisions(:) ! model decisions ! define input-output variables - integer(i4b),intent(inout) :: iFile ! index of current forcing file in forcing file list - integer(i4b),intent(inout) :: iRead ! index of read position in time dimension in current netcdf file - integer(i4b),intent(inout) :: ncid ! netcdf file identifier + integer(i4b),intent(inout) :: iFile ! index of current forcing file in forcing file list + integer(i4b),intent(inout) :: iRead ! index of read position in time dimension in current netcdf file + integer(i4b),intent(inout) :: ncid ! netcdf file identifier ! define output variables - integer(i4b),intent(out) :: time_data(:) ! vector of time data for a given time step - type(gru_hru_double) :: forcStruct ! x%gru(:)%hru(:)%var(:) -- model forcing data - integer(i4b),intent(out) :: err ! error code - character(*),intent(out) :: message ! error message + integer(i4b),intent(out) :: time_data(:) ! vector of time data for a given time step + type(gru_hru_double) :: forcStruct ! x%gru(:)%hru(:)%var(:) -- model forcing data + integer(i4b),intent(out) :: err ! error code + character(*),intent(out) :: message ! error message ! define local variables - integer(i4b) :: nHRUlocal ! number of HRUs in the local simulation - integer(i4b) :: iGRU,iHRU ! index of GRU and HRU - character(len=256),save :: infile ! filename - character(len=256) :: cmessage ! error message for downwind routine - real(rkind) :: startJulDay ! julian day at the start of the year - real(rkind) :: currentJulDay ! Julian day of current time step + integer(i4b) :: nRemain ! number of steps remaining in the simulation + integer(i4b) :: nData ! number of steps remaining in the file + integer(i4b) :: nRead ! number of steps in the data read + integer(i4b) :: nHRUlocal ! number of HRUs in the local simulation + integer(i4b) :: iline ! loop through lines in the file + integer(i4b) :: iTime ! index of time + integer(i4b) :: jRead ! index of time in data subset + integer(i4b) :: iGRU,iHRU ! index of GRU and HRU + character(len=256),save :: infile ! filename + real(rkind) :: dsec ! double precision seconds (not used) + real(rkind) :: dataJulDay ! julian day of current forcing data step being read + real(rkind) :: startJulDay ! julian day at the start of the year + real(rkind) :: currentJulDay ! Julian day of current time step + logical(lgt) :: isNewFile ! .true. if reading a new forcing file + logical(lgt) :: isRead ! .true. if reading data logical(lgt),parameter :: checkTime=.false. ! flag to check the time + ! error control + integer(i4b) :: ierr ! local error code + character(len=256) :: cmessage ! error message for downwind routine ! Start procedure here err=0; message="read_force/" + ! initialize new file + isNewFile = .false. + ! get the number of HRUs in the local simulation nHRUlocal = sum(gru_struc(:)%hruCount) @@ -131,12 +159,14 @@ subroutine read_force(istep,iFile,iRead,ncid,time_data,forcStruct,err,message) call getFirstTimestep(currentJulDay,iFile,iRead,ncid,err,cmessage) if(err/=0)then; message=trim(message)//trim(cmessage); return; end if + ! flag new file + isNewFile=.true. + end if ! if the file is not yet open ! ********************************************************************************************** - ! ***** part 1: if file open, check to see if we've reached the end of the file, if so close it, - ! ***** and open new file - ! ***** Then read the data + ! ***** part 1a: if file open, check to see if we've reached the end of the file, if so close it, + ! ***** and open new file ! ********************************************************************************************** if(ncid>0)then @@ -161,31 +191,114 @@ subroutine read_force(istep,iFile,iRead,ncid,time_data,forcStruct,err,message) call openForcingFile(iFile,trim(infile),ncId,err,cmessage) if(err/=0)then; message=trim(message)//trim(cmessage); return; end if + ! flag new file + isNewFile=.true. + ! reset iRead since we opened a new file iRead=1 end if ! if we've passed the end of the NetCDF file - ! read forcing data - call readForcingData(currentJulDay,ncId,iFile,iRead,nHRUlocal,time_data,forcStruct,err,message) - if(err/=0)then; message=trim(message)//trim(cmessage); return; end if - ! check that the file was in fact open else message=trim(message)//'expect the file to be open' err=20; return end if ! end ncid open check + ! ********************************************************************************************** + ! ***** part 1b: allocate space for data structures and read the data + ! ********************************************************************************************** + + ! re-define data length and arrays if new file + if(isNewFile)then + + ! deallocate space + if(allocated(fulltimeVec)) deallocate(fulltimeVec) + if(allocated(fullforcingStruct)) deallocate(fullforcingStruct) + + ! get the number of time steps to read + select case(model_decisions(iLookDECISIONS%read_force)%iDecision) + case(readPerStep); nRead=1 ! ** read forcing data per time step (defualt) + case(readFullSeries) ! ** read full forcing series + nRemain = (numtim - iStep) + 1 ! number of remaining time steps in simulation + nData = (forcFileInfo(iFile)%nTimeSteps - iRead) + 1 ! number of remaining data steps in file (starting with iRead) + nRead = min(nRemain,nData) ! number of data steps to read + case default; err=10; message=trim(message)//'unable to identify option to read forcing data'; return + end select + + ! allocate space + allocate(fulltimeVec(nRead), stat=ierr) + if(ierr/=0)then; err=20; message=trim(message)//'problem allocating space for fulltimeVec'; return; endif + allocate(fullforcingStruct(nRead), source=forcStruct, stat=ierr) + if(ierr/=0)then; err=20; message=trim(message)//'problem allocating space for fullforcingStruct'; return; endif + + ! define starting index for the data read + ixStartRead = iRead + + endif ! if new file + + ! update ixStartRead and nRead for readPerStep + if(model_decisions(iLookDECISIONS%read_force)%iDecision == readPerStep)then + ixStartRead = iRead + nRead = 1 + endif + + ! check if we need to read the data + isRead = (isNewFile .or. model_decisions(iLookDECISIONS%read_force)%iDecision == readPerStep) + + ! read forcing data + ! NOTE: reads data into global variables fulltimeVec and fullforcingStruct + if(isRead)then + call readForcingData(ncId,iFile,ixStartRead,nRead,nHRUlocal,err,cmessage) + if(err/=0)then; message=trim(message)//trim(cmessage); return; end if + endif ! of reading the forcing data + + ! get the index in the data structures + jRead = (iRead - ixStartRead) + 1 + + ! get forcing structure for the desired time + forcStruct = fullforcingStruct(jRead) + !print*, forcStruct%gru(1)%hru(1)%var + #endif ! ********************************************************************************************** ! ***** part 2: compute time ! ********************************************************************************************** + ! check that the computed julian day matches the time information in the NetCDF file + dataJulDay = fulltimeVec(jRead)/forcFileInfo(iFile)%convTime2Days + refJulDay_data + if(abs(currentJulDay - dataJulDay) > timeDiffTol)then + write(message,'(a,f18.8,a,f18.8)') trim(message)//'date for time step: ',dataJulDay,' differs from the expected date: ',currentJulDay + err=40; return + end if + + ! convert julian day to time vector + ! NOTE: use small offset to force ih=0 at the start of the day + call compcalday(dataJulDay+smallOffset, & ! input = julian day + time_data(iLookTIME%iyyy), & ! output = year + time_data(iLookTIME%im), & ! output = month + time_data(iLookTIME%id), & ! output = day + time_data(iLookTIME%ih), & ! output = hour + time_data(iLookTIME%imin),dsec, & ! output = minute/second + err,cmessage) ! output = error control + if(err/=0)then; message=trim(message)//trim(cmessage); return; end if + + ! check to see if any of the time data is missing -- note that it is OK if ih_tz or imin_tz are missing + if((time_data(iLookTIME%iyyy)==integerMissing) .or. & + (time_data(iLookTIME%im) ==integerMissing) .or. & + (time_data(iLookTIME%id) ==integerMissing) .or. & + (time_data(iLookTIME%ih) ==integerMissing) .or. & + (time_data(iLookTIME%imin)==integerMissing) )then + do iline=1,size(time_data) + if(time_data(iline)==integerMissing)then; err=40; message=trim(message)//"variableMissing[var='"//trim(time_meta(iline)%varname)//"']"; return; end if + end do + end if ! if time data is missing + ! compute the julian day at the start of the year - call compjulday(time_data(iLookTIME%iyyy), & ! input = year - 1, 1, 1, 1, 0._rkind, & ! input = month, day, hour, minute, second - startJulDay,err,cmessage) ! output = julian day (fraction of day) + error control + call compjulday(time_data(iLookTIME%iyyy), & ! input = year + 1, 1, 1, 1, 0._rkind, & ! input = month, day, hour, minute, second + startJulDay,err,cmessage) ! output = julian day (fraction of day) + error control if(err/=0)then; message=trim(message)//trim(cmessage); return; end if ! compute the fractional julian day for the current time step @@ -414,74 +527,39 @@ end subroutine openForcingFile ! ************************************************************************* ! * read the NetCDF forcing data ! ************************************************************************* - subroutine readForcingData(currentJulDay,ncId,iFile,iRead,nHRUlocal,time_data,forcStruct,err,message) + subroutine readForcingData(ncId,iFile,ixStartRead,nRead,nHRUlocal,err,message) USE netcdf ! netcdf capability USE time_utils_module,only:compcalday ! convert julian day to calendar date USE time_utils_module,only:compJulDay ! convert calendar date to julian day USE get_ixname_module,only:get_ixForce ! identify index of named variable ! dummy variables - real(rkind),intent(in) :: currentJulDay ! Julian day of current time step - integer(i4b) ,intent(in) :: ncId ! NetCDF ID - integer(i4b) ,intent(in) :: iFile ! index of forcing file - integer(i4b) ,intent(in) :: iRead ! index in data file - integer(i4b) ,intent(in) :: nHRUlocal ! number of HRUs in the local simulation - integer(i4b),intent(out) :: time_data(:) ! vector of time data for a given time step - type(gru_hru_double) :: forcStruct ! x%gru(:)%hru(:)%var(:) -- model forcing data - integer(i4b) ,intent(out) :: err ! error code - character(*) ,intent(out) :: message ! error message + integer(i4b) ,intent(in) :: ncId ! NetCDF ID + integer(i4b) ,intent(in) :: iFile ! index of forcing file + integer(i4b) ,intent(in) :: ixStartRead ! starting index in data file + integer(i4b) ,intent(in) :: nRead ! number of time steps for the local data read + integer(i4b) ,intent(in) :: nHRUlocal ! number of HRUs in the local simulation + integer(i4b) ,intent(out) :: err ! error code + character(*) ,intent(out) :: message ! error message ! local variables - character(len=256) :: cmessage ! error message for downwind routine - integer(i4b) :: varId ! variable identifier - character(len = nf90_max_name) :: varName ! dimenison name - real(rkind) :: varTime(1) ! time variable of current forcing data step being read + integer(i4b) :: varId ! variable identifier + character(len = nf90_max_name) :: varName ! dimenison name ! other local variables - integer(i4b) :: iGRU,iHRU ! index of GRU and HRU - integer(i4b) :: iHRU_global ! index of HRU in the NetCDF file - integer(i4b) :: iHRU_local ! index of HRU in the data subset - integer(i4b) :: iline ! loop through lines in the file - integer(i4b) :: iNC ! loop through variables in forcing file - integer(i4b) :: iVar ! index of forcing variable in forcing data vector - real(rkind) :: dsec ! double precision seconds (not used) - real(rkind) :: dataJulDay ! julian day of current forcing data step being read - real(rkind),dimension(nHRUlocal) :: dataVec ! vector of data - real(rkind),dimension(1) :: dataVal ! single data value - real(rkind),parameter :: dataMin=-1._rkind ! minimum allowable data value (all forcing variables should be positive) - logical(lgt),dimension(size(forc_meta)) :: checkForce ! flags to check forcing data variables exist - logical(lgt),parameter :: simultaneousRead=.true. ! flag to denote reading all HRUs at once + integer(i4b) :: iTime ! time index + integer(i4b) :: iGRU,iHRU ! index of GRU and HRU + integer(i4b) :: iHRU_global ! index of HRU in the NetCDF file + integer(i4b) :: iHRU_local ! index of HRU in the data subset + integer(i4b) :: iline ! loop through lines in the file + integer(i4b) :: iNC ! loop through variables in forcing file + integer(i4b) :: iVar ! index of forcing variable in forcing data vector + real(rkind),dimension(nHRUlocal,nRead) :: dataMatrix ! vector of data + real(rkind),parameter :: dataMin=-1._rkind ! minimum allowable data value (all forcing variables should be positive) + logical(lgt),dimension(size(forc_meta)) :: checkForce ! flags to check forcing data variables exist ! Start procedure here err=0; message="readForcingData/" - ! initialize time and forcing data structures - time_data(:) = integerMissing - - ! read time data from iRead location in netcdf file - err = nf90_inq_varid(ncid,'time',varId); if(err/=nf90_noerr)then; message=trim(message)//'trouble finding time variable/'//trim(nf90_strerror(err)); return; endif - err = nf90_get_var(ncid,varId,varTime,start=(/iRead/)); if(err/=nf90_noerr)then; message=trim(message)//'trouble reading time variable/'//trim(nf90_strerror(err)); return; endif - - ! check that the computed julian day matches the time information in the NetCDF file - dataJulDay = varTime(1)/forcFileInfo(iFile)%convTime2Days + refJulDay_data - if(abs(currentJulDay - dataJulDay) > timeDiffTol)then - write(message,'(a,f18.8,a,f18.8)') trim(message)//'date for time step: ',dataJulDay,' differs from the expected date: ',currentJulDay - err=40; return - end if - - ! convert julian day to time vector - ! NOTE: use small offset to force ih=0 at the start of the day - call compcalday(dataJulDay+smallOffset, & ! input = julian day - time_data(iLookTIME%iyyy), & ! output = year - time_data(iLookTIME%im), & ! output = month - time_data(iLookTIME%id), & ! output = day - time_data(iLookTIME%ih), & ! output = hour - time_data(iLookTIME%imin),dsec, & ! output = minute/second - err,cmessage) ! output = error control - if(err/=0)then; message=trim(message)//trim(cmessage); return; end if - - ! check to see if any of the time data is missing -- note that it is OK if ih_tz or imin_tz are missing - if((time_data(iLookTIME%iyyy)==integerMissing) .or. (time_data(iLookTIME%im)==integerMissing) .or. (time_data(iLookTIME%id)==integerMissing) .or. (time_data(iLookTIME%ih)==integerMissing) .or. (time_data(iLookTIME%imin)==integerMissing))then - do iline=1,size(time_data) - if(time_data(iline)==integerMissing)then; err=40; message=trim(message)//"variableMissing[var='"//trim(time_meta(iline)%varname)//"']"; return; end if - end do - end if + ! read time data starting from iRead location in netcdf file + err = nf90_inq_varid(ncid,'time',varId); if(err/=nf90_noerr)then; message=trim(message)//'trouble finding time variable/'//trim(nf90_strerror(err)); return; endif + err = nf90_get_var(ncid,varId,fulltimeVec,start=(/ixStartRead/),count=(/nRead/)); if(err/=nf90_noerr)then; message=trim(message)//'trouble reading time variable/'//trim(nf90_strerror(err)); return; endif ! initialize flags for forcing data checkForce(:) = .false. @@ -501,46 +579,40 @@ subroutine readForcingData(currentJulDay,ncId,iFile,iRead,nHRUlocal,time_data,fo err=nf90_inquire_variable(ncid,iNC,name=varName) if(err/=nf90_noerr)then; message=trim(message)//'problem reading forcing variable name from netCDF: '//trim(nf90_strerror(err)); return; endif - ! read forcing data for all HRUs - if(simultaneousRead)then - err=nf90_get_var(ncid,forcFileInfo(iFile)%data_id(ivar),dataVec,start=(/ixHRUfile_min,iRead/),count=(/nHRUlocal,1/)) - if(err/=nf90_noerr)then; message=trim(message)//'problem reading forcing data: '//trim(varName)//'/'//trim(nf90_strerror(err)); return; endif - endif - - ! loop through GRUs and HRUs - do iGRU=1,size(gru_struc) - do iHRU=1,gru_struc(iGRU)%hruCount - - ! define global HRU - iHRU_global = gru_struc(iGRU)%hruInfo(iHRU)%hru_nc - iHRU_local = (iHRU_global - ixHRUfile_min)+1 - - ! read forcing data for a single HRU - if(.not.simultaneousRead)then - err=nf90_get_var(ncid,forcFileInfo(iFile)%data_id(ivar),dataVal,start=(/iHRU_global,iRead/)) - if(err/=nf90_noerr)then; message=trim(message)//'problem reading forcing data: '//trim(varName)//'/'//trim(nf90_strerror(err)); return; endif - endif - - ! check the number of HRUs - if(iHRU_global > nHRUfile)then - message=trim(message)//'HRU index exceeds the number of HRUs in the forcing data file' - err=20; return - endif - - ! get individual data value - if(simultaneousRead) dataVal(1) = dataVec(iHRU_local) - - ! check individual data value - if(dataVal(1) nHRUfile)then + message=trim(message)//'HRU index exceeds the number of HRUs in the forcing data file' + err=20; return + endif + + ! check individual data value + if(dataMatrix(iHRU_local,iTime) < dataMin)then + write(message,'(a,f13.5)') trim(message)//'forcing data for variable '//trim(varname)//' is less than minimum allowable value ', dataMin + err=20; return + endif + + ! put the data into structures + fullforcingStruct(iTime)%gru(iGRU)%hru(iHRU)%var(ivar) = dataMatrix(iHRU_local,iTime) + + end do ! looping through HRUs within a given GRU + end do ! looping through GRUs + + end do ! looping through time end do ! loop through forcing variables diff --git a/build/source/netcdf/def_output.f90 b/build/source/netcdf/def_output.f90 index 217debb26..0153bf30a 100644 --- a/build/source/netcdf/def_output.f90 +++ b/build/source/netcdf/def_output.f90 @@ -24,6 +24,7 @@ module def_output_module USE netcdf_util_module,only:nc_file_close ! close NetCDF files USE f2008funcs_module,only:cloneStruc ! used to "clone" data structures -- temporary replacement of the intrinsic allocate(a, source=b) USE nrtype, integerMissing=>nr_integerMissing ! top-level data types +USE globalData, only: numtim USE globalData, only: outputPrecision ! data structure for output precision USE globalData, only: chunkSize ! size of chunks to write USE globalData, only: outputCompressionLevel ! netcdf deflate level @@ -314,11 +315,14 @@ subroutine def_variab(ncid,iFreq,spatialDesire,timeDesire,metaData,ivtype,err,me if (metaData(iVar)%statIndex(iFreq)==integerMissing .and. metaData(iVar)%varName/='time') cycle ! ---------- get the dimension IDs (use cloneStruc, given source) ---------- - gruChunk = min(nGRUrun, chunkSize) - hruChunk = min(nHRUrun, chunkSize) + gruChunk = min(nGRUrun, chunkSize) + hruChunk = min(nHRUrun, chunkSize) timeChunk = chunkSize layerChunk = 1 + ! set the chunk size to the number of time steps for single basin runs + if(nGRUrun == 1) timeChunk = numtim + ! special case of the time variable if(metaData(iVar)%varName == 'time')then call cloneStruc(dimensionIDs, lowerBound=1, source=(/Timestep_DimID/),err=err,message=cmessage); writechunk=(/ timeChunk /) diff --git a/build/source/netcdf/modelwrite.f90 b/build/source/netcdf/modelwrite.f90 index b59f9fb3a..b6fd537e9 100644 --- a/build/source/netcdf/modelwrite.f90 +++ b/build/source/netcdf/modelwrite.f90 @@ -31,6 +31,9 @@ module modelwrite_module USE globalData,only: integerMissing, realMissing ! provide access to global data +USE globalData,only:nGRUrun ! number of GRUs in the run +USE globalData,only:nHRUrun ! number of HRUs in the run +USE globalData,only:maxLayers ! maximum number of layers USE globalData,only:gru_struc ! gru->hru mapping structure ! provide access to the derived types to define the data structures @@ -68,7 +71,6 @@ module modelwrite_module private public::writeParm public::writeData -public::writeBasin public::writeTime public::writeRestart ! define dimension lengths @@ -142,26 +144,26 @@ subroutine writeParm(ispatial,struct,meta,err,message) end subroutine writeParm ! ************************************************************************************** - ! public subroutine writeData: write model time-dependent data + ! public subroutine writeData: write model time-dependent data for each HRU ! ************************************************************************************** - subroutine writeData(finalizeStats,outputTimestep,nHRUrun,maxLayers,meta,stat,dat,map,indx,err,message) + subroutine writeData(is_bufferedWrite,finalizeStats,outputTimestep,maxWrite,stat,timestepData,meta,map,indx,err,message) USE data_types,only:var_info ! metadata type USE var_lookup,only:maxVarStat ! index into stats structure USE var_lookup,only:iLookVarType ! index into type structure USE var_lookup,only:iLookINDEX ! index into index structure - USE var_lookup,only:iLookSTAT ! index into stat structure + USE var_lookup,only:iLookFREQ ! index into freq structure USE globalData,only:outFreq,ncid ! output file information 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 implicit none ! declare dummy variables + logical(lgt) ,intent(in) :: is_bufferedWrite ! flag for buffered write logical(lgt) ,intent(in) :: finalizeStats(:) ! flags to finalize statistics integer(i4b) ,intent(in) :: outputTimestep(:) ! output time step - integer(i4b) ,intent(in) :: nHRUrun ! number of HRUs in the run domain - integer(i4b) ,intent(in) :: maxLayers ! maximum number of layers + integer(i4b) ,intent(in) :: maxWrite ! maximum number of steps written type(var_info),intent(in) :: meta(:) ! meta data class(*) ,intent(in) :: stat ! stats data - class(*) ,intent(in) :: dat ! timestep data + class(*) ,intent(in) :: timestepData(:) ! timestep data integer(i4b) ,intent(in) :: map(:) ! map into stats child struct type(gru_hru_intVec) ,intent(in) :: indx ! index data integer(i4b) ,intent(out) :: err ! error code @@ -172,15 +174,19 @@ subroutine writeData(finalizeStats,outputTimestep,nHRUrun,maxLayers,meta,stat,da integer(i4b) :: iVar ! variable index integer(i4b) :: iStat ! statistics index integer(i4b) :: iFreq ! frequency index + integer(i4b) :: iTime ! time index integer(i4b) :: ncVarID ! used only for time integer(i4b) :: nSnow ! number of snow layers integer(i4b) :: nSoil ! number of soil layers integer(i4b) :: nLayers ! total number of layers + integer(i4b) :: ixStart ! index of the start of data write + integer(i4b) :: nSpace ! number of spatial data elements ! output arrays integer(i4b) :: datLength ! length of each data vector integer(i4b) :: maxLength ! maximum length of each data vector - real(rkind) :: realVec(nHRUrun) ! real vector for all HRUs in the run domain - real(rkind) :: realArray(nHRUrun,maxLayers+1) ! real array for all HRUs in the run domain + real(rkind) :: timeBuffer(maxWrite) ! buffer for all time steps + real(rkind) :: realBuffer(nHRUrun,maxWrite) ! buffer for all HRUs in the run domain + time steps + real(rkind) :: realArray(nHRUrun,maxLayers+1) ! real array for all HRUs in the run domain integer(i4b) :: intArray(nHRUrun,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 @@ -194,28 +200,61 @@ subroutine writeData(finalizeStats,outputTimestep,nHRUrun,maxLayers,meta,stat,da ! skip frequencies that are not needed if(.not.outFreq(iFreq)) cycle + ! restrict attention to the timestep data if buffered write + if(is_bufferedWrite .and. iFreq/=iLookFREQ%timestep) cycle + ! check that we have finalized statistics for a given frequency if(.not.finalizeStats(iFreq)) cycle + ! get the start index + if(is_bufferedWrite)then + ixStart = 1 + else + ixStart = outputTimestep(iFreq) + endif + ! loop through model variables do iVar = 1,size(meta) + !if(meta(iVar)%varDesire) print*, meta(iVar)%varName + + ! **************************************************************************** + ! *** write time information -- instantaneous + ! **************************************************************************** + ! handle time first if (meta(iVar)%varName=='time')then + ! get variable index err = nf90_inq_varid(ncid(iFreq),trim(meta(iVar)%varName),ncVarID) call netcdf_err(err,message); if (err/=0) return + ! define HRUs and GRUs (only write once) iGRU=1; iHRU=1 - ! data bound write - select type(dat) ! forcStruc + + ! data bound array access + select type(timestepData) ! forcStruc class is (gru_hru_double) ! x%gru(:)%hru(:)%var(:) - err = nf90_put_var(ncid(iFreq),ncVarID,(/dat%gru(iGRU)%hru(iHRU)%var(iVar)/),start=(/outputTimestep(iFreq)/),count=(/1/)) - call netcdf_err(err,message); if (err/=0) return - cycle ! move onto the next variable + + ! put data in time buffer + do iTime=1,maxWrite + timeBuffer(iTime) = timestepData(iTime)%gru(iGRU)%hru(iHRU)%var(iVar) + end do + + ! check we found the class class default; err=20; message=trim(message)//'time variable must be of type gru_hru_double (forcing data structure)'; return - end select - end if ! id time + end select ! type of data structure + + ! write time + err = nf90_put_var(ncid(iFreq),ncVarID,(/timeBuffer/),start=(/ixStart/),count=(/maxWrite/)) + call netcdf_err(err,message); if (err/=0) return + cycle ! move onto the next variable + + end if ! if time + + ! **************************************************************************** + ! *** write scalar variables + ! **************************************************************************** ! define the statistics index iStat = meta(iVar)%statIndex(iFreq) @@ -225,30 +264,109 @@ subroutine writeData(finalizeStats,outputTimestep,nHRUrun,maxLayers,meta,stat,da ! stats output: only scalar variable type if(meta(iVar)%varType==iLookVarType%scalarv) then - select type(stat) - class is (gru_hru_doubleVec) - ! loop through HRUs and GRUs, and place data in the single vector + ! ----- writing buffered output data --------------------------------------- + if(is_bufferedWrite)then + + ! loop through time, HRUs and GRUs, and place data in the buffer + do iTime=1,maxWrite do iGRU=1,size(gru_struc) + + ! identify data structures + select type(timestepData) + + ! *** HRU structures (indices) + class is (gru_hru_int) + if(iGRU==1) nSpace = nHRUrun + do iHRU=1,gru_struc(iGRU)%hruCount + realBuffer(gru_struc(iGRU)%hruInfo(iHRU)%hru_ix,iTime) = timestepData(iTime)%gru(iGRU)%hru(iHRU)%var(map(iVar)) + end do ! hru + + ! *** HRU structures (forcing, prognostic, diagnostic, ...) + class is (gru_hru_double) + if(iGRU==1) nSpace = nHRUrun + do iHRU=1,gru_struc(iGRU)%hruCount + realBuffer(gru_struc(iGRU)%hruInfo(iHRU)%hru_ix,iTime) = timestepData(iTime)%gru(iGRU)%hru(iHRU)%var(map(iVar)) + end do ! hru + + ! *** GRU structures (basin-average variables, ...) + class is (gru_double) + if(iGRU==1) nSpace = nGRUrun + realBuffer(iGRU,iTime) = timestepData(iTime)%gru(iGRU)%var(map(iVar)) + + class default; err=20; message=trim(message)//'cannot identify data class'; return + end select + + end do ! gru + end do ! time + + ! write data -- note that the number of GRUs is less than or equal to the number of HRUs, so pass using nSpace + err = nf90_put_var(ncid(iFreq),meta(iVar)%ncVarID(iFreq),realBuffer(1:nSpace,1:maxWrite),start=(/1,1/),count=(/nSpace,maxWrite/)) + call netcdf_err(err,message); if (err/=0) return + + ! ----- writing statistics ------------------------------------------------- + + ! check that we are not writing buffered data -- writing statistics + else + + ! check that maxWrite==1 + if(maxWrite/=1)then + message=trim(message)//'expect maxWrite=1 when not writing buffered output' + err=20; return + endif + + ! loop through GRUs + do iGRU=1,size(gru_struc) + + ! identify data structures + select type(stat) + + ! *** HRU structures (forcing, prognostic, diagnostic, ...) + class is (gru_hru_doubleVec) + if(iGRU==1) nSpace = nHRUrun do iHRU=1,gru_struc(iGRU)%hruCount - realVec(gru_struc(iGRU)%hruInfo(iHRU)%hru_ix) = stat%gru(iGRU)%hru(iHRU)%var(map(iVar))%dat(iFreq) + realBuffer(gru_struc(iGRU)%hruInfo(iHRU)%hru_ix,1) = stat%gru(iGRU)%hru(iHRU)%var(map(iVar))%dat(iFreq) end do - end do + + ! *** GRU structures (basin-average variables, ...) + class is (gru_doubleVec) + if(iGRU==1) nSpace = nGRUrun + realBuffer(iGRU,1) = stat%gru(iGRU)%var(map(iVar))%dat(iFreq) + + ! check statistics type + class default + message=trim(message)//'stats must be scalarv and of type gru_hru_doubleVec or gru_doubleVec' + err=20; return + end select ! stat + + end do ! (looping through GRUs - ! write data - err = nf90_put_var(ncid(iFreq),meta(iVar)%ncVarID(iFreq),realVec,start=(/1,outputTimestep(iFreq)/),count=(/nHRUrun,1/)) + ! write data + err = nf90_put_var(ncid(iFreq),meta(iVar)%ncVarID(iFreq),realBuffer(1:nSpace,1),start=(/1,outputTimestep(iFreq)/),count=(/nSpace,1/)) + call netcdf_err(err,message); if (err/=0) return - class default; err=20; message=trim(message)//'stats must be scalarv and of type gru_hru_doubleVec'; return - end select ! stat + endif ! (if not buffered write -- statistics) + + ! **************************************************************************** + ! *** write non-scalar variables (regular data structures -- instantaneous) + ! **************************************************************************** ! non-scalar variables: regular data structures else + ! cannot get here if buffered write + if(is_bufferedWrite)then + message=trim(message)//'cannot write non-scalar variables in buffered write ['//trim(meta(iVar)%varName)//']' + err=20; return + endif + ! initialize the data vectors - select type (dat) + select type (timestepData) class is (gru_hru_doubleVec); realArray(:,:) = realMissing; dataType=ixReal class is (gru_hru_intVec); intArray(:,:) = integerMissing; dataType=ixInteger - class default; err=20; message=trim(message)//'data must not be scalarv and either of type gru_hru_doubleVec or gru_hru_intVec'; return + class default + message=trim(message)//'data must not be scalarv and either of type gru_hru_doubleVec or gru_hru_intVec ['//trim(meta(iVar)%varName)//']' + err=20; return end select ! loop thru GRUs and HRUs @@ -273,9 +391,9 @@ subroutine writeData(finalizeStats,outputTimestep,nHRUrun,maxLayers,meta,stat,da end select ! vartype ! get the data vectors - select type (dat) - class is (gru_hru_doubleVec); realArray(gru_struc(iGRU)%hruInfo(iHRU)%hru_ix,1:datLength) = dat%gru(iGRU)%hru(iHRU)%var(iVar)%dat(:) - class is (gru_hru_intVec); intArray(gru_struc(iGRU)%hruInfo(iHRU)%hru_ix,1:datLength) = dat%gru(iGRU)%hru(iHRU)%var(iVar)%dat(:) + select type (timestepData) + class is (gru_hru_doubleVec); realArray(gru_struc(iGRU)%hruInfo(iHRU)%hru_ix,1:datLength) = timestepData(1)%gru(iGRU)%hru(iHRU)%var(iVar)%dat(:) + class is (gru_hru_intVec); intArray(gru_struc(iGRU)%hruInfo(iHRU)%hru_ix,1:datLength) = timestepData(1)%gru(iGRU)%hru(iHRU)%var(iVar)%dat(:) 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 @@ -312,77 +430,6 @@ subroutine writeData(finalizeStats,outputTimestep,nHRUrun,maxLayers,meta,stat,da end subroutine writeData - ! ************************************************************************************** - ! public subroutine writeBasin: write basin-average variables - ! ************************************************************************************** - subroutine writeBasin(iGRU,finalizeStats,outputTimestep,meta,stat,dat,map,err,message) - USE data_types,only:var_info ! metadata type - USE var_lookup,only:maxVarStat ! index into stats structure - USE var_lookup,only:iLookVarType ! index into type structure - USE globalData,only:outFreq,ncid ! output file information - 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 - implicit none - - ! declare dummy variables - integer(i4b) ,intent(in) :: iGRU ! GRU index - logical(lgt) ,intent(in) :: finalizeStats(:) ! flags to finalize statistics - integer(i4b) ,intent(in) :: outputTimestep(:) ! output time step - type(var_info),intent(in) :: meta(:) ! meta data - type(dlength) ,intent(in) :: stat(:) ! stats data - type(dlength) ,intent(in) :: dat(:) ! timestep data - integer(i4b) ,intent(in) :: map(:) ! map into stats child struct - integer(i4b) ,intent(out) :: err ! error code - character(*) ,intent(out) :: message ! error message - ! local variables - integer(i4b) :: iVar ! variable index - integer(i4b) :: iStat ! statistics index - integer(i4b) :: iFreq ! frequency index - ! initialize error control - err=0;message="f-writeBasin/" - - ! loop through output frequencies - do iFreq=1,maxvarFreq - - ! skip frequencies that are not needed - if(.not.outFreq(iFreq)) cycle - - ! check that we have finalized statistics for a given frequency - if(.not.finalizeStats(iFreq)) cycle - - ! loop through model variables - do iVar = 1,size(meta) - - ! define the statistics index - iStat = meta(iVar)%statIndex(iFreq) - - ! check that the variable is desired - if (iStat==integerMissing.or.trim(meta(iVar)%varName)=='unknown') cycle - - ! stats/data output - select data type - select case (meta(iVar)%varType) - - case (iLookVarType%scalarv) - err = nf90_put_var(ncid(iFreq),meta(iVar)%ncVarID(iFreq),(/stat(map(iVar))%dat(iFreq)/),start=(/iGRU,outputTimestep(iFreq)/),count=(/1,1/)) - - case (iLookVarType%routing) - if (iFreq==1 .and. outputTimestep(iFreq)==1) then - err = nf90_put_var(ncid(iFreq),meta(iVar)%ncVarID(iFreq),(/dat(iVar)%dat/),start=(/1/),count=(/1000/)) - end if - - case default - err=40; message=trim(message)//"unknownVariableType[name='"//trim(meta(iVar)%varName)//"';type='"//trim(get_varTypeName(meta(iVar)%varType))// "']"; return - end select ! variable type - - ! process error code - if (err.ne.0) message=trim(message)//trim(meta(iVar)%varName)//'_'//trim(get_statName(iStat)) - call netcdf_err(err,message); if (err/=0) return - - end do ! iVar - end do ! iFreq - - end subroutine writeBasin - ! ************************************************************************************** ! public subroutine writeTime: write current time to all files ! **************************************************************************************