From 415fe898a7d496959950dbd18c87a3505146820f Mon Sep 17 00:00:00 2001 From: Kyle Klenk <kyle.c.klenk@gmail.com> Date: Thu, 25 Aug 2022 21:40:09 +0000 Subject: [PATCH] implemented new subroutine read_param_all_hru --- .../file_access_actor_subroutine_wrappers.hpp | 2 + build/makefile | 3 +- .../file_access_actor/file_access_actor.cpp | 12 +- .../file_access_actor/read_param_all_hru.f90 | 275 +++++++++ build/source/driver/SummaActors_setup.f90 | 13 +- build/source/driver/summaActors_util.f90 | 226 ------- build/source/engine/read_paramActors.f90 | 572 +++++++++--------- 7 files changed, 568 insertions(+), 535 deletions(-) create mode 100644 build/source/actors/file_access_actor/read_param_all_hru.f90 diff --git a/build/includes/file_access_actor/file_access_actor_subroutine_wrappers.hpp b/build/includes/file_access_actor/file_access_actor_subroutine_wrappers.hpp index ade26ec..a510b92 100644 --- a/build/includes/file_access_actor/file_access_actor_subroutine_wrappers.hpp +++ b/build/includes/file_access_actor/file_access_actor_subroutine_wrappers.hpp @@ -33,4 +33,6 @@ extern "C" { void Write_HRU_Param(void* handle_ncid, int* indxGRU, int* indxHRU, int* err); void readAttributeFileAccessActor(int* num_gru, int* err); + + void readParamFileAccessActor(int* start_gru, int* num_gru, int* err); } diff --git a/build/makefile b/build/makefile index 3d6283a..e2c5668 100644 --- a/build/makefile +++ b/build/makefile @@ -144,7 +144,8 @@ SUMMA_FILEACCESS_INTERFACE = \ initOutputStruc.f90 \ deallocateOutputStruc.f90 \ cppwrap_fileAccess.f90 \ - read_attribute_all_hru.f90 + read_attribute_all_hru.f90 \ + read_param_all_hru.f90 FILEACCESS_INTERFACE = $(patsubst %, $(FILE_ACCESS_DIR)/%, $(SUMMA_FILEACCESS_INTERFACE)) diff --git a/build/source/actors/file_access_actor/file_access_actor.cpp b/build/source/actors/file_access_actor/file_access_actor.cpp index 560dfdf..7086276 100644 --- a/build/source/actors/file_access_actor/file_access_actor.cpp +++ b/build/source/actors/file_access_actor/file_access_actor.cpp @@ -240,7 +240,7 @@ void initalizeFileAccessActor(stateful_actor<file_access_state>* self) { Init_OutputStruct(self->state.handle_forcing_file_info, &self->state.outputStrucSize, &self->state.numGRU, &self->state.err); - // Read In all of the attribres for the number of GRUs in the run Domian + // Read In all of the attribres for the number of GRUs in the run Domian readAttributeFileAccessActor(&self->state.numGRU, &err); if (err != 0) { aout(self) << "ERROR: FILE_ACCESS_ACTOR readAttributeFilAccessActor() \n"; @@ -250,6 +250,16 @@ void initalizeFileAccessActor(stateful_actor<file_access_state>* self) { return; } + // Read in all of the parmeters for the number of GRUs in the run Domain + readParamFileAccessActor(&self->state.startGRU, &self->state.numGRU, &err); + if (err != 0) { + aout(self) << "ERROR: FILE_ACCESS_ACTOR readParamFileAccessActor() \n"; + std::string function = "readParamFileAccessActor"; + self->send(self->state.parent, file_access_actor_err_v, function); + self->quit(); + return; + } + // Initalize the output manager self->state.output_manager = new OutputManager(self->state.num_vectors_in_output_manager, self->state.numGRU); diff --git a/build/source/actors/file_access_actor/read_param_all_hru.f90 b/build/source/actors/file_access_actor/read_param_all_hru.f90 new file mode 100644 index 0000000..31ff21b --- /dev/null +++ b/build/source/actors/file_access_actor/read_param_all_hru.f90 @@ -0,0 +1,275 @@ +module read_param_all_hru + USE, intrinsic :: iso_c_binding + USE nrtype + implicit none + private + public::read_param_file_access_actor +contains +subroutine read_param_file_access_actor(startGRU,num_gru,err) bind(C, name="readParamFileAccessActor") + ! used to read model initial conditions + USE summaActors_FileManager,only:SETTINGS_PATH ! path for metadata files + USE summaActors_FileManager,only:PARAMETER_TRIAL ! file with parameter trial values + USE get_ixname_module,only:get_ixparam,get_ixbpar ! access function to find index of elements in structure + USE globalData,only:index_map,gru_struc ! mapping from global HRUs to the elements in the data structures + USE var_lookup,only:iLookPARAM,iLookTYPE,iLookID ! named variables to index elements of the data vectors + USE globalData,only:integerMissing ! missing integer + USE globalData,only:realMissing ! missing real number + + USE netcdf + USE netcdf_util_module,only:nc_file_close ! close netcdf file + USE netcdf_util_module,only:nc_file_open ! open netcdf file + USE netcdf_util_module,only:netcdf_err ! netcdf error handling function + + USE globalData,only:outputStructure + + implicit none + ! define input + + integer(c_int),intent(in) :: startGRU ! starting Index of gru Batch + integer(c_int),intent(in) :: num_gru ! number of GRUs in the run_domain + integer(c_int),intent(out) :: err ! error code + + ! define local variables + character(len=256) :: message ! error message + character(len=1024) :: cmessage ! error message for downwind routine + character(LEN=1024) :: infile ! input filename + integer(i4b) :: localHRU_ix ! index of HRU within data structure + integer(i4b) :: ixParam ! index of the model parameter in the data structure + ! indices/metadata in the NetCDF file + integer(i4b) :: ncid ! netcdf id + integer(i4b) :: nDims ! number of dimensions + integer(i4b) :: nVars ! number of variables + integer(i4b) :: idimid ! dimension index + integer(i4b) :: ivarid ! variable index + character(LEN=64) :: dimName ! dimension name + character(LEN=64) :: parName ! parameter name + integer(i4b) :: dimLength ! dimension length + integer(i4b) :: nHRU_file ! number of HRUs in the parafile + integer(i4b) :: nGRU_file ! number of GRUs in the parafile + integer(i4b) :: nSoil_file ! number of soil layers in the file + integer(i4b) :: idim_list(2) ! list of dimension ids + ! data in the netcdf file + integer(i4b) :: parLength ! length of the parameter data + integer(8),allocatable :: hruId(:) ! HRU identifier in the file + real(dp),allocatable :: parVector(:) ! model parameter vector + logical :: fexist ! inquire whether the parmTrial file exists + integer(i4b) :: fHRU ! index of HRU in input file + integer(i4b) :: iGRU + integer(i4b) :: iHRU + + err=0; message="read_param_all_hru.f90/" + + + ! ********************************************************************************************** + ! * open files, etc. + ! ********************************************************************************************** + + infile = trim(SETTINGS_PATH)//trim(PARAMETER_TRIAL) ! build filename + + ! check whether the user-specified file exists and warn if it does not + inquire(file=trim(infile),exist=fexist) + if (.not.fexist) then + write(*,'(A)') NEW_LINE('A')//'!! WARNING: trial parameter file not found; proceeding instead with other default parameters; check path in file manager input if this was not the desired behavior'//NEW_LINE('A') + return + endif + + ! open trial parameters file if it exists + call nc_file_open(trim(infile),nf90_nowrite,ncid,err,cmessage) + if(err/=0)then; message=trim(message)//trim(cmessage); return; end if + + ! get the number of variables in the parameter file + err=nf90_inquire(ncid, nDimensions=nDims, nVariables=nVars) + call netcdf_err(err,message); if (err/=0) then; err=20; return; end if + + ! initialize the number of HRUs + nHRU_file=integerMissing + nGRU_file=integerMissing + + ! get the length of the dimensions + do idimid=1,nDims + ! get the dimension name and length + err=nf90_inquire_dimension(ncid, idimid, name=dimName, len=dimLength) + if(err/=0)then; message=trim(message)//trim(cmessage); return; end if + ! get the number of HRUs + if(trim(dimName)=='hru') nHRU_file=dimLength + if(trim(dimName)=='gru') nGRU_file=dimLength + end do + + ! allocate hruID vector + allocate(hruId(nHRU_file)) + + ! check HRU dimension exists + if(nHRU_file==integerMissing)then + message=trim(message)//'unable to identify HRU dimension in file '//trim(infile) + err=20; return + endif + ! ********************************************************************************************** + ! * read the HRU index + ! ********************************************************************************************** + ! loop through the parameters in the NetCDF file + do ivarid=1,nVars + + ! get the parameter name + err=nf90_inquire_variable(ncid, ivarid, name=parName) + call netcdf_err(err,message) + if (err/=0) then + err=20 + print*, message + return + end if + + ! special case of the HRU id + if(trim(parName)=='hruIndex' .or. trim(parName)=='hruId')then + + ! read HRUs + err=nf90_get_var(ncid, ivarid, hruId) + if(err/=0)then + message=trim(message)//trim(cmessage) + print*, message + return + end if + + endif ! if the HRU id + + end do ! looping through variables in the file + + ! ********************************************************************************************** + ! * read the local parameters and the basin parameters + ! ********************************************************************************************** + do ivarid=1,nVars + ! get the parameter name + err=nf90_inquire_variable(ncid, ivarid, name=parName) + call netcdf_err(err,message); if (err/=0) then; err=20; return; end if + + ! get the local parameters + ixParam = get_ixparam( trim(parName) ) + if(ixParam/=integerMissing)then + ! ********************************************************************************************** + ! * read the local parameters + ! ********************************************************************************************** + + ! get the variable shape + err=nf90_inquire_variable(ncid, ivarid, nDims=nDims, dimids=idim_list) + if(err/=0)then + message=trim(message)//trim(cmessage) + print*, message + return + end if + + ! get the length of the depth dimension (if it exists) + if(nDims==2)then + ! get the information on the 2nd dimension for 2-d variables + err=nf90_inquire_dimension(ncid, idim_list(2), dimName, nSoil_file) + if(err/=0)then + message=trim(message)//trim(cmessage) + print*, message + return + end if + + ! check that it is the depth dimension + if(trim(dimName)/='depth')then + message=trim(message)//'expect 2nd dimension of 2-d variable to be depth (dimension name = '//trim(dimName)//')' + err=20; return + endif + + ! TODO: implement this check + ! ! check that the dimension length is correct + ! if(size(outputStructur(1)%mparStruct%var(ixParam)%dat) /= nSoil_file)then + ! message=trim(message)//'unexpected number of soil layers in parameter file' + ! err=20; return + ! endif + + ! define parameter length + parLength = nSoil_file + else + parLength = 1 + endif ! if two dimensions + + ! allocate space for model parameters + allocate(parVector(parLength),stat=err) + if(err/=0)then + message=trim(message)//'problem allocating space for parameter vector' + err=20; return + endif + + ! loop through all hrus + do iGRU=1, num_gru + do iHRU=1, gru_struc(iGRU)%hruCount + localHRU_ix=index_map(iHRU)%localHRU_ix + fHRU = gru_struc(iGRU)%hruInfo(localHRU_ix)%hru_nc + ! read parameter data + select case(nDims) + case(1); err=nf90_get_var(ncid, ivarid, parVector, start=(/fHRU/), count=(/1/) ) + case(2); err=nf90_get_var(ncid, ivarid, parVector, start=(/fHRU,1/), count=(/1,nSoil_file/) ) + case default; err=20; message=trim(message)//'unexpected number of dimensions for parameter '//trim(parName) + end select + + ! error check for the parameter read + if(err/=0)then + message=trim(message)//trim(cmessage) + print*, message + return + end if + + ! populate parameter structures + select case(nDims) + case(1); outputStructure(1)%mparStruct(1)%gru(iGRU)%hru(localHRU_ix)%var(ixParam)%dat(:) = parVector(1) ! also distributes scalar across depth dimension + case(2); outputStructure(1)%mparStruct(1)%gru(iGRU)%hru(localHRU_ix)%var(ixParam)%dat(:) = parVector(:) + case default; err=20; message=trim(message)//'unexpected number of dimensions for parameter '//trim(parName) + end select + end do + end do + + ! deallocate space for model parameters + deallocate(parVector,stat=err) + if(err/=0)then + message=trim(message)//'problem deallocating space for parameter vector' + print*, message + err=20; return + endif + + ! ********************************************************************************************** + ! * read the basin parameters + ! ********************************************************************************************** + + ! get the basin parameters + else + ! get the parameter index + ixParam = get_ixbpar( trim(parName) ) + + ! allow extra variables in the file that are not used + if(ixParam==integerMissing) cycle + + ! allocate space for model parameters + allocate(parVector(nGRU_file),stat=err) + if(err/=0)then + message=trim(message)//'problem allocating space for parameter vector' + print*, message + err=20; return + endif + + ! read parameter data + err=nf90_get_var(ncid, ivarid, parVector ) + if(err/=0)then; message=trim(message)//trim(cmessage); return; end if + + ! populate parameter structures + do iGRU=1, num_gru + outputStructure(1)%bparStruct(1)%gru(1)%var(ixParam) = parVector(iGRU+startGRU-1) + end do + + ! deallocate space for model parameters + deallocate(parVector,stat=err) + if(err/=0)then + message=trim(message)//'problem deallocating space for parameter vector' + print*, message + err=20; return + endif + endif + + end do ! (looping through the parameters in the NetCDF file) + + ! Now we must close the netcdf file + call nc_file_close(ncid,err,message) + if(err/=0)then;message=trim(message)//trim(cmessage);return;end if +end subroutine read_param_file_access_actor +end module read_param_all_hru \ No newline at end of file diff --git a/build/source/driver/SummaActors_setup.f90 b/build/source/driver/SummaActors_setup.f90 index 8b268e5..7789608 100755 --- a/build/source/driver/SummaActors_setup.f90 +++ b/build/source/driver/SummaActors_setup.f90 @@ -94,7 +94,7 @@ contains USE read_attribute_module,only:read_attribute ! module to read local attributes USE paramCheck_module,only:paramCheck ! module to check consistency of model parameters USE pOverwrite_module,only:pOverwrite ! module to overwrite default parameter values with info from the Noah tables - USE read_param4chm_module,only:read_param4chm ! module to read model parameter sets + USE read_param4chm_module,only:read_param ! module to read model parameter sets USE ConvE2Temp_module,only:E2T_lookup ! module to calculate a look-up table for the temperature-enthalpy conversion USE var_derive_module,only:fracFuture ! module to calculate the fraction of runoff in future time steps (time delay histogram) USE module_sf_noahmplsm,only:read_mp_veg_parameters ! module to read NOAH vegetation tables @@ -188,21 +188,12 @@ contains ! ***************************************************************************** ! *** read local attributes for each HRU ! ***************************************************************************** - - - ! read local attributes for each HRU call read_attribute(indxHRU,indxGRU,attrStruct,typeStruct,idStruct,err,cmessage) if(err/=0)then; message=trim(message)//trim(cmessage); return; endif - ! ***************************************************************************** - ! *** read Noah vegetation and soil tables - ! ***************************************************************************** - ! define monthly fraction of green vegetation greenVegFrac_monthly = (/0.01_dp, 0.02_dp, 0.03_dp, 0.07_dp, 0.50_dp, 0.90_dp, 0.95_dp, 0.96_dp, 0.65_dp, 0.24_dp, 0.11_dp, 0.02_dp/) - - ! define urban vegetation category select case(trim(model_decisions(iLookDECISIONS%vegeParTbl)%cDecision)) case('USGS'); urbanVegCategory = 1 @@ -236,7 +227,7 @@ contains ! ***************************************************************************** ! *** read trial model parameter values for each HRU, and populate initial data structures ! ***************************************************************************** - call read_param4chm(indxHRU,indxGRU,iRunMode,startGRU, & + call read_param(indxHRU,indxGRU,iRunMode,startGRU, & mparStruct,bparStruct,err,cmessage) if(err/=0)then; message=trim(message)//trim(cmessage); return; endif ! ***************************************************************************** diff --git a/build/source/driver/summaActors_util.f90 b/build/source/driver/summaActors_util.f90 index 5150565..0b76107 100755 --- a/build/source/driver/summaActors_util.f90 +++ b/build/source/driver/summaActors_util.f90 @@ -36,236 +36,10 @@ implicit none private ! routines to make public -! public::getCommandArguments4chm public::stop_program public::handle_err contains - ! ************************************************************************************************** - ! * obtain the command line arguments - ! ************************************************************************************************** -! subroutine getCommandArguments4chm(file_manager_path,summaFileManagerFile,err,message) -! ! data types -! USE summa4chm_type, only:summa4chm_type_dec ! master summa data type -! ! provide access to named parameters -! USE globalData,only:iRunModeFull,iRunModeGRU,iRunModeHRU -! USE globalData,only:ixProgress_it,ixProgress_im,ixProgress_id,ixProgress_ih,ixProgress_never -! USE globalData,only:ixRestart_iy,ixRestart_im,ixRestart_id,ixRestart_end,ixRestart_never -! USE globalData,only:noNewFiles,newFileEveryOct1 -! ! provide access to runtime options -! USE globalData,only: startGRU ! index of the starting GRU for parallelization run -! USE globalData,only: checkHRU ! index of the HRU for a single HRU run -! USE globalData,only: iRunMode ! define the current running mode -! USE globalData,only: newOutputFile ! define option for new output file -! USE globalData,only: ixProgress ! define frequency to write progress -! USE globalData,only: ixRestart ! define frequency to write restart files -! USE globalData,only: output_fileSuffix ! suffix for the output file -! implicit none -! ! dummy variables -! character(len=256),intent(in) :: file_manager_path ! path/name of file defining directories and files -! character(len=256),intent(inout) :: summaFileManagerFile ! path/name of file defining directories and files -! integer(i4b),intent(out) :: err ! error code -! character(*),intent(out) :: message ! error message -! ! local variables -! integer(i4b) :: nGRU ! number of grouped response units -! integer(i4b) :: nHRU ! number of global hydrologic response units -! integer(i4b) :: iArgument ! index of command line argument -! integer(i4b) :: nArgument ! number of command line arguments -! character(len=256),allocatable :: argString(:) ! string to store command line arguments -! integer(i4b) :: nLocalArgument ! number of command line arguments to read for a switch -! character(len=70), parameter :: spaces = '' ! setting a blank string -! ! version information generated during compiling -! INCLUDE 'summaversion.inc' -! ! --------------------------------------------------------------------------------------- -! ! initialize error control -! err=0; message='getCommandArguments4chm/' - -! nArgument = 6 -! allocate(argString(nArgument)) -! argString(1) = '-p' -! argString(2) = 'never' -! argString(3) = '-s' -! argString(4) = '_testSumma' -! argString(5) = '-m' -! argString(6) = file_manager_path - -! ! initialize command line argument variables -! startGRU = integerMissing; checkHRU = integerMissing -! nGRU = integerMissing; nHRU = integerMissing -! newOutputFile = noNewFiles -! iRunMode = iRunModeFull - -! ! loop through all command arguments -! nLocalArgument = 0 -! do iArgument = 1,nArgument -! if (nLocalArgument>0) then; nLocalArgument = nLocalArgument -1; cycle; end if ! skip the arguments have been read -! select case (trim(argString(iArgument))) - -! case ('-m', '--master') -! ! update arguments -! nLocalArgument = 1 -! if (iArgument+nLocalArgument>nArgument)then -! message="missing argument file_suffix; type 'summa.exe --help' for correct usage" -! err=1; return -! endif -! ! get name of master control file -! summaFileManagerFile=trim(argString(iArgument+1)) -! print "(A)", "file_master is '"//trim(summaFileManagerFile)//"'." - -! ! define the formation of new output files -! case ('-n', '--newFile') -! ! check that the number of command line arguments is correct -! nLocalArgument = 1 ! expect just one argument for new output files -! if (iArgument+nLocalArgument>nArgument)then -! message="missing argument file_suffix; type 'summa.exe --help' for correct usage" -! err=1; return -! endif -! ! get the decision for the formation of new output files -! select case( trim(argString(iArgument+1)) ) -! case('noNewFiles'); newOutputFile = noNewFiles -! case('newFileEveryOct1'); newOutputFile = newFileEveryOct1 -! case default -! message='unknown option for new output file: expect "noNewFiles" or "newFileEveryOct1"' -! err=1; return -! end select - -! case ('-s', '--suffix') -! ! define file suffix -! nLocalArgument = 1 -! ! check if the number of command line arguments is correct -! if (iArgument+nLocalArgument>nArgument) then -! message="missing argument file_suffix; type 'summa.exe --help' for correct usage" -! err=1; return -! endif -! output_fileSuffix=trim(argString(iArgument+1)) -! print "(A)", "file_suffix is '"//trim(output_fileSuffix)//"'." - -! case ('-h', '--hru') -! ! define a single HRU run -! if (iRunMode == iRunModeGRU)then -! message="single-HRU run and GRU-parallelization run cannot be both selected." -! err=1; return -! endif -! iRunMode=iRunModeHRU -! nLocalArgument = 1 -! ! check if the number of command line arguments is correct -! if (iArgument+nLocalArgument>nArgument) call handle_err(1,"missing argument checkHRU; type 'summa.exe --help' for correct usage") -! read(argString(iArgument+1),*) checkHRU ! read the index of the HRU for a single HRU run -! nHRU=1; nGRU=1 ! nHRU and nGRU are both one in this case -! ! examines the checkHRU is correct -! if (checkHRU<1) then -! message="illegal iHRU specification; type 'summa.exe --help' for correct usage" -! err=1; return -! else -! print '(A)',' Single-HRU run activated. HRU '//trim(argString(iArgument+1))//' is selected for simulation.' -! end if - -! case ('-g','--gru') -! ! define a GRU parallelization run; get the starting GRU and countGRU -! if (iRunMode == iRunModeHRU)then -! message="single-HRU run and GRU-parallelization run cannot be both selected." -! err=1; return -! endif -! iRunMode=iRunModeGRU -! nLocalArgument = 2 -! ! check if the number of command line arguments is correct -! if (iArgument+nLocalArgument>nArgument)then -! message="missing argument startGRU or countGRU; type 'summa.exe --help' for correct usage" -! err=1; return -! endif -! read(argString(iArgument+1),*) startGRU ! read the argument of startGRU -! read(argString(iArgument+2),*) nGRU ! read the argument of countGRU -! if (startGRU<1 .or. nGRU<1) then -! message='startGRU and countGRU must be larger than 1.' -! err=1; return -! else -! print '(A)', ' GRU-Parallelization run activated. '//trim(argString(iArgument+2))//' GRUs are selected for simulation.' -! end if - -! case ('-p', '--progress') -! ! define the frequency to print progress -! nLocalArgument = 1 -! ! check if the number of command line arguments is correct -! if (iArgument+nLocalArgument>nArgument)then -! message="missing argument freqProgress; type 'summa.exe --help' for correct usage" -! err=1; return -! endif -! select case (trim(argString(iArgument+1))) -! case ('t' , 'timestep'); ixProgress = ixProgress_it -! case ('h' , 'hour'); ixProgress = ixProgress_ih -! case ('d' , 'day'); ixProgress = ixProgress_id ! default -! case ('m' , 'month'); ixProgress = ixProgress_im -! case ('n' , 'never'); ixProgress = ixProgress_never -! case default -! message='unknown frequency to print progress' -! err=1; return -! end select - -! case ('-r', '--restart') -! ! define the frequency to write restart files -! nLocalArgument = 1 -! ! check if the number of command line arguments is correct -! if (iArgument+nLocalArgument>nArgument)then -! message="missing argument freqRestart; type 'summa.exe --help' for correct usage" -! err=1; return -! endif -! select case (trim(argString(iArgument+1))) -! case ('y' , 'year'); ixRestart = ixRestart_iy -! case ('m' , 'month'); ixRestart = ixRestart_im -! case ('d' , 'day'); ixRestart = ixRestart_id -! case ('e' , 'end'); ixRestart = ixRestart_end -! case ('n' , 'never'); ixRestart = ixRestart_never -! case default -! message='unknown frequency to write restart files' -! err=1; return -! end select - -! ! do nothing -! case ('-v','--version') - -! ! print help message -! case ('--help') -! call printCommandHelp - -! case default -! call printCommandHelp -! message='unknown command line option' -! err=1; return - -! end select -! end do ! looping through command line arguments - -! ! check if master_file has been received. -! if (len(trim(summaFileManagerFile))==0)then -! message="master_file is not received; type 'summa.exe --help' for correct usage" -! err=1; return -! endif - -! ! set startGRU for full run -! if (iRunMode==iRunModeFull) startGRU=1 - -! end subroutine getCommandArguments4chm - - ! ************************************************************************************************** - ! print the correct command line usage of SUMMA - ! ************************************************************************************************** -! subroutine printCommandHelp() -! implicit none -! ! command line usage -! print "(//A)",'Usage: summa.exe -m master_file [-s fileSuffix] [-g startGRU countGRU] [-h iHRU] [-r freqRestart] [-p freqProgress] [-c]' -! print "(A,/)", ' summa.exe summa executable' -! print "(A)", 'Running options:' -! print "(A)", ' -m --master Define path/name of master file (required)' -! print "(A)", ' -n --newFile Define frequency [noNewFiles,newFileEveryOct1] of new output files' -! print "(A)", ' -s --suffix Add fileSuffix to the output files' -! print "(A)", ' -g --gru Run a subset of countGRU GRUs starting from index startGRU' -! print "(A)", ' -h --hru Run a single HRU with index of iHRU' -! print "(A)", ' -r --restart Define frequency [y,m,d,e,never] to write restart files' -! print "(A)", ' -p --progress Define frequency [m,d,h,never] to print progress' -! print "(A)", ' -v --version Display version information of the current build' -! stop -! end subroutine printCommandHelp - ! ************************************************************************************************** ! error handler ! ************************************************************************************************** diff --git a/build/source/engine/read_paramActors.f90 b/build/source/engine/read_paramActors.f90 index bac913f..93f4a95 100755 --- a/build/source/engine/read_paramActors.f90 +++ b/build/source/engine/read_paramActors.f90 @@ -41,319 +41,299 @@ USE data_types,only:var_dlength ! spatial double data type: x%gru(:)%hru(: implicit none private -public::read_param4chm +public::read_param contains - ! ************************************************************************************************ - ! public subroutine read_param4chm: read trial model parameter values - ! ************************************************************************************************ - subroutine read_param4chm(indxHRU,indxGRU,iRunMode,startGRU,mparStruct,bparStruct,err,message) - ! used to read model initial conditions - USE summaActors_FileManager,only:SETTINGS_PATH ! path for metadata files - USE summaActors_FileManager,only:PARAMETER_TRIAL ! file with parameter trial values - USE get_ixname_module,only:get_ixparam,get_ixbpar ! access function to find index of elements in structure - USE globalData,only:index_map,gru_struc ! mapping from global HRUs to the elements in the data structures - USE var_lookup,only:iLookPARAM,iLookTYPE,iLookID ! named variables to index elements of the data vectors - implicit none - ! define input - integer(i4b),intent(in) :: indxHRU - integer(i4b),intent(in) :: indxGRU - integer(i4b),intent(in) :: iRunMode ! run mode - integer(i4b),intent(in) :: startGRU ! index of single GRU if runMode = startGRU -! type(var_i8),intent(in) :: idStruct ! local labels for hru and gru IDs - ! define output - type(var_dlength),intent(inout) :: mparStruct ! model parameters - type(var_d),intent(inout) :: bparStruct ! basin parameters - integer(i4b),intent(out) :: err ! error code - character(*),intent(out) :: message ! error message - ! define local variables - character(len=1024) :: cmessage ! error message for downwind routine - character(LEN=1024) :: infile ! input filename - integer(i4b) :: localHRU_ix ! index of HRU within data structure - integer(i4b) :: ixParam ! index of the model parameter in the data structure - ! indices/metadata in the NetCDF file - integer(i4b) :: ncid ! netcdf id - integer(i4b) :: nDims ! number of dimensions - integer(i4b) :: nVars ! number of variables - integer(i4b) :: idimid ! dimension index - integer(i4b) :: ivarid ! variable index - character(LEN=64) :: dimName ! dimension name - character(LEN=64) :: parName ! parameter name - integer(i4b) :: dimLength ! dimension length - integer(i4b) :: nHRU_file ! number of HRUs in the parafile - integer(i4b) :: nGRU_file ! number of GRUs in the parafile - integer(i4b) :: nSoil_file ! number of soil layers in the file - integer(i4b) :: idim_list(2) ! list of dimension ids - ! data in the netcdf file - integer(i4b) :: parLength ! length of the parameter data - integer(8),allocatable :: hruId(:) ! HRU identifier in the file - real(dp),allocatable :: parVector(:) ! model parameter vector - logical :: fexist ! inquire whether the parmTrial file exists - integer(i4b) :: fHRU ! index of HRU in input file - - ! Start procedure here - err=0; message="read_param4chm/" - - ! ********************************************************************************************** - ! * open files, etc. - ! ********************************************************************************************** - - ! build filename - infile = trim(SETTINGS_PATH)//trim(PARAMETER_TRIAL) - - ! check whether the user-specified file exists and warn if it does not - inquire(file=trim(infile),exist=fexist) - if (.not.fexist) then - write(*,'(A)') NEW_LINE('A')//'!! WARNING: trial parameter file not found; proceeding instead with other default parameters; check path in file manager input if this was not the desired behavior'//NEW_LINE('A') - return - endif - - ! open trial parameters file if it exists - call nc_file_open(trim(infile),nf90_nowrite,ncid,err,cmessage) - if(err/=0)then; message=trim(message)//trim(cmessage); return; end if - - ! get the number of variables in the parameter file - err=nf90_inquire(ncid, nDimensions=nDims, nVariables=nVars) - call netcdf_err(err,message); if (err/=0) then; err=20; return; end if - - ! initialize the number of HRUs - nHRU_file=integerMissing - nGRU_file=integerMissing - - ! get the length of the dimensions - do idimid=1,nDims - ! get the dimension name and length - err=nf90_inquire_dimension(ncid, idimid, name=dimName, len=dimLength) - if(err/=0)then; message=trim(message)//trim(cmessage); return; end if - ! get the number of HRUs - if(trim(dimName)=='hru') nHRU_file=dimLength - if(trim(dimName)=='gru') nGRU_file=dimLength - end do - - ! allocate hruID vector - allocate(hruId(nHRU_file)) - - ! check HRU dimension exists - if(nHRU_file==integerMissing)then - message=trim(message)//'unable to identify HRU dimension in file '//trim(infile) - err=20; return - endif - -! check have the correct number of HRUs -! if ((irunMode==irunModeFull).and.(nHRU_file/=nHRU)) then -! message=trim(message)//'incorrect number of HRUs in file '//trim(infile) -! err=20; return -! endif -! if ((irunMode==irunModeHRU).and.(nHRU_file<checkHRU)) then -! message=trim(message)//'not enough HRUs in file '//trim(infile) -! err=20; return -! endif - -! ! check have the correct number of GRUs -! if ((irunMode==irunModeGRU).and.(nGRU_file<startGRU).and.(nGRU_file/=integerMissing)) then -! message=trim(message)//'not enough GRUs in file '//trim(infile) -! err=20; return -! endif -! if ((irunMode==irunModeFull).and.(nGRU_file/=nGRU).and.(nGRU_file/=integerMissing)) then -! message=trim(message)//'incorrect number of GRUs in file '//trim(infile) -! err=20; return -! endif - - ! ********************************************************************************************** - ! * read the HRU index - ! ********************************************************************************************** - - ! loop through the parameters in the NetCDF file - do ivarid=1,nVars - - ! get the parameter name - err=nf90_inquire_variable(ncid, ivarid, name=parName) - call netcdf_err(err,message); if (err/=0) then; err=20; return; end if - - ! special case of the HRU id - if(trim(parName)=='hruIndex' .or. trim(parName)=='hruId')then - - ! read HRUs - err=nf90_get_var(ncid, ivarid, hruId) - if(err/=0)then; message=trim(message)//trim(cmessage); return; end if - - ! check HRUs -- expect HRUs to be in the same order as the local attributes -! if (iRunMode==iRunModeFull) then -! !iGRU=index_map(indxHRU)%gru_ix -! localHRU_ix=index_map(indxHRU)%localHRU_ix -! if((hruId(indxHRU)>0).and.(hruId(indxHRU)/=idStruct%var(iLookID%hruId)))then -! write(message,'(a,i0,a,i0,a)') trim(message)//'mismatch for HRU ', idStruct%var(iLookID%hruId), '(param HRU = ', hruId(indxHRU), ')' -! err=20; return -! endif - -! else if (iRunMode==iRunModeGRU) then -! ! do iHRU=1,nHRU -! !iGRU=index_map(indxHRU)%gru_ix -! localHRU_ix=index_map(indxHRU)%localHRU_ix -! fHRU = gru_struc(indxGRU)%hruInfo(localHRU_ix)%hru_nc -! if(hruId(fHRU)/=idStruct%var(iLookID%hruId))then -! write(message,'(a,i0,a,i0,a)') trim(message)//'mismatch for HRU ', idStruct%var(iLookID%hruId), '(param HRU = ', hruId(indxHRU), ')' -! err=20; return -! endif -! ! enddo - -! else if (iRunMode==iRunModeHRU) then -! !iGRU=index_map(1)%gru_ix -! localHRU_ix=index_map(indxHRU)%localHRU_ix -! if(hruId(checkHRU)/=idStruct%var(iLookID%hruId))then -! write(message,'(a,i0,a,i0,a)') trim(message)//'mismatch for HRU ', idStruct%var(iLookID%hruId), '(param HRU = ', hruId(indxHRU), ')' -! err=20; return -! endif - - ! error check -! else -! err = 20; message = 'run mode not recognized'; return; -! end if - - endif ! if the HRU id - - end do ! looping through variables in the file - - ! ********************************************************************************************** - ! * read the local parameters and the basin parameters - ! ********************************************************************************************** - - ! loop through the parameters in the NetCDF file - do ivarid=1,nVars - - ! get the parameter name - err=nf90_inquire_variable(ncid, ivarid, name=parName) - call netcdf_err(err,message); if (err/=0) then; err=20; return; end if - - ! get the local parameters - ixParam = get_ixparam( trim(parName) ) - if(ixParam/=integerMissing)then +! ************************************************************************************************ +! public subroutine read_param4chm: read trial model parameter values +! ************************************************************************************************ +subroutine read_param(indxHRU,indxGRU,iRunMode,startGRU,mparStruct,bparStruct,err,message) + ! used to read model initial conditions + USE summaActors_FileManager,only:SETTINGS_PATH ! path for metadata files + USE summaActors_FileManager,only:PARAMETER_TRIAL ! file with parameter trial values + USE get_ixname_module,only:get_ixparam,get_ixbpar ! access function to find index of elements in structure + USE globalData,only:index_map,gru_struc ! mapping from global HRUs to the elements in the data structures + USE var_lookup,only:iLookPARAM,iLookTYPE,iLookID ! named variables to index elements of the data vectors + implicit none + ! define input + integer(i4b),intent(in) :: indxHRU + integer(i4b),intent(in) :: indxGRU + integer(i4b),intent(in) :: iRunMode ! run mode + integer(i4b),intent(in) :: startGRU ! index of single GRU if runMode = startGRU + ! type(var_i8),intent(in) :: idStruct ! local labels for hru and gru IDs + ! define output + type(var_dlength),intent(inout) :: mparStruct ! model parameters + type(var_d),intent(inout) :: bparStruct ! basin parameters + integer(i4b),intent(out) :: err ! error code + character(*),intent(out) :: message ! error message + ! define local variables + character(len=1024) :: cmessage ! error message for downwind routine + character(LEN=1024) :: infile ! input filename + integer(i4b) :: localHRU_ix ! index of HRU within data structure + integer(i4b) :: ixParam ! index of the model parameter in the data structure + ! indices/metadata in the NetCDF file + integer(i4b) :: ncid ! netcdf id + integer(i4b) :: nDims ! number of dimensions + integer(i4b) :: nVars ! number of variables + integer(i4b) :: idimid ! dimension index + integer(i4b) :: ivarid ! variable index + character(LEN=64) :: dimName ! dimension name + character(LEN=64) :: parName ! parameter name + integer(i4b) :: dimLength ! dimension length + integer(i4b) :: nHRU_file ! number of HRUs in the parafile + integer(i4b) :: nGRU_file ! number of GRUs in the parafile + integer(i4b) :: nSoil_file ! number of soil layers in the file + integer(i4b) :: idim_list(2) ! list of dimension ids + ! data in the netcdf file + integer(i4b) :: parLength ! length of the parameter data + integer(8),allocatable :: hruId(:) ! HRU identifier in the file + real(dp),allocatable :: parVector(:) ! model parameter vector + logical :: fexist ! inquire whether the parmTrial file exists + integer(i4b) :: fHRU ! index of HRU in input file + + ! Start procedure here + err=0; message="read_param4chm/" ! ********************************************************************************************** - ! * read the local parameters + ! * open files, etc. ! ********************************************************************************************** - ! get the variable shape - err=nf90_inquire_variable(ncid, ivarid, nDims=nDims, dimids=idim_list) - if(err/=0)then; message=trim(message)//trim(cmessage); return; end if - - ! get the length of the depth dimension (if it exists) - if(nDims==2)then - - ! get the information on the 2nd dimension for 2-d variables - err=nf90_inquire_dimension(ncid, idim_list(2), dimName, nSoil_file) - if(err/=0)then; message=trim(message)//trim(cmessage); return; end if - - ! check that it is the depth dimension - if(trim(dimName)/='depth')then - message=trim(message)//'expect 2nd dimension of 2-d variable to be depth (dimension name = '//trim(dimName)//')' - err=20; return - endif - - ! check that the dimension length is correct - if(size(mparStruct%var(ixParam)%dat) /= nSoil_file)then - message=trim(message)//'unexpected number of soil layers in parameter file' - err=20; return - endif - - ! define parameter length - parLength = nSoil_file - - else - parLength = 1 - endif ! if two dimensions - - ! allocate space for model parameters - allocate(parVector(parLength),stat=err) - if(err/=0)then - message=trim(message)//'problem allocating space for parameter vector' - err=20; return - endif - - ! map to the GRUs and HRUs - !iGRU=index_map(indxHRU)%gru_ix - localHRU_ix=index_map(indxHRU)%localHRU_ix - fHRU = gru_struc(indxGRU)%hruInfo(localHRU_ix)%hru_nc - - ! read parameter data - select case(nDims) - case(1); err=nf90_get_var(ncid, ivarid, parVector, start=(/fHRU/), count=(/1/) ) - case(2); err=nf90_get_var(ncid, ivarid, parVector, start=(/fHRU,1/), count=(/1,nSoil_file/) ) - case default; err=20; message=trim(message)//'unexpected number of dimensions for parameter '//trim(parName) - end select - - ! error check for the parameter read - if(err/=0)then; message=trim(message)//trim(cmessage); return; end if - - ! populate parameter structures - select case(nDims) - case(1); mparStruct%var(ixParam)%dat(:) = parVector(1) ! also distributes scalar across depth dimension - case(2); mparStruct%var(ixParam)%dat(:) = parVector(:) - case default; err=20; message=trim(message)//'unexpected number of dimensions for parameter '//trim(parName) - end select - - ! end do ! looping through HRUs - - ! deallocate space for model parameters - deallocate(parVector,stat=err) - if(err/=0)then - message=trim(message)//'problem deallocating space for parameter vector' - err=20; return - endif - - ! ********************************************************************************************** - ! * read the basin parameters - ! ********************************************************************************************** - - ! get the basin parameters - else - - ! get the parameter index - ixParam = get_ixbpar( trim(parName) ) - - ! allow extra variables in the file that are not used - if(ixParam==integerMissing) cycle + ! build filename + infile = trim(SETTINGS_PATH)//trim(PARAMETER_TRIAL) - ! allocate space for model parameters - allocate(parVector(nGRU_file),stat=err) - if(err/=0)then - message=trim(message)//'problem allocating space for parameter vector' - err=20; return + ! check whether the user-specified file exists and warn if it does not + inquire(file=trim(infile),exist=fexist) + if (.not.fexist) then + write(*,'(A)') NEW_LINE('A')//'!! WARNING: trial parameter file not found; proceeding instead with other default parameters; check path in file manager input if this was not the desired behavior'//NEW_LINE('A') + return endif - ! read parameter data - err=nf90_get_var(ncid, ivarid, parVector ) + ! open trial parameters file if it exists + call nc_file_open(trim(infile),nf90_nowrite,ncid,err,cmessage) if(err/=0)then; message=trim(message)//trim(cmessage); return; end if - ! populate parameter structures - if (iRunMode==iRunModeGRU) then - !do iGRU=1,nGRU - bparStruct%var(ixParam) = parVector(indxGRU+startGRU-1) - !end do ! looping through GRUs - else if (iRunMode==iRunModeFull) then - !do iGRU=1,nGRU - bparStruct%var(ixParam) = parVector(indxGRU) - !end do ! looping through GRUs - else if (iRunMode==iRunModeHRU) then - err = 20; message='checkHRU run mode not working'; return; + ! get the number of variables in the parameter file + err=nf90_inquire(ncid, nDimensions=nDims, nVariables=nVars) + call netcdf_err(err,message); if (err/=0) then; err=20; return; end if + + ! initialize the number of HRUs + nHRU_file=integerMissing + nGRU_file=integerMissing + + ! get the length of the dimensions + do idimid=1,nDims + ! get the dimension name and length + err=nf90_inquire_dimension(ncid, idimid, name=dimName, len=dimLength) + if(err/=0)then; message=trim(message)//trim(cmessage); return; end if + ! get the number of HRUs + if(trim(dimName)=='hru') nHRU_file=dimLength + if(trim(dimName)=='gru') nGRU_file=dimLength + end do + + ! allocate hruID vector + allocate(hruId(nHRU_file)) + + ! check HRU dimension exists + if(nHRU_file==integerMissing)then + message=trim(message)//'unable to identify HRU dimension in file '//trim(infile) + err=20; return endif + - ! deallocate space for model parameters - deallocate(parVector,stat=err) - if(err/=0)then - message=trim(message)//'problem deallocating space for parameter vector' - err=20; return - endif + ! ********************************************************************************************** + ! * read the HRU index + ! ********************************************************************************************** + ! loop through the parameters in the NetCDF file + do ivarid=1,nVars + + ! get the parameter name + err=nf90_inquire_variable(ncid, ivarid, name=parName) + call netcdf_err(err,message); if (err/=0) then; err=20; return; end if + + ! special case of the HRU id + if(trim(parName)=='hruIndex' .or. trim(parName)=='hruId')then + + ! read HRUs + err=nf90_get_var(ncid, ivarid, hruId) + if(err/=0)then; message=trim(message)//trim(cmessage); return; end if + + ! check HRUs -- expect HRUs to be in the same order as the local attributes + ! if (iRunMode==iRunModeFull) then + ! !iGRU=index_map(indxHRU)%gru_ix + ! localHRU_ix=index_map(indxHRU)%localHRU_ix + ! if((hruId(indxHRU)>0).and.(hruId(indxHRU)/=idStruct%var(iLookID%hruId)))then + ! write(message,'(a,i0,a,i0,a)') trim(message)//'mismatch for HRU ', idStruct%var(iLookID%hruId), '(param HRU = ', hruId(indxHRU), ')' + ! err=20; return + ! endif + + ! else if (iRunMode==iRunModeGRU) then + ! ! do iHRU=1,nHRU + ! !iGRU=index_map(indxHRU)%gru_ix + ! localHRU_ix=index_map(indxHRU)%localHRU_ix + ! fHRU = gru_struc(indxGRU)%hruInfo(localHRU_ix)%hru_nc + ! if(hruId(fHRU)/=idStruct%var(iLookID%hruId))then + ! write(message,'(a,i0,a,i0,a)') trim(message)//'mismatch for HRU ', idStruct%var(iLookID%hruId), '(param HRU = ', hruId(indxHRU), ')' + ! err=20; return + ! endif + ! ! enddo + + ! else if (iRunMode==iRunModeHRU) then + ! !iGRU=index_map(1)%gru_ix + ! localHRU_ix=index_map(indxHRU)%localHRU_ix + ! if(hruId(checkHRU)/=idStruct%var(iLookID%hruId))then + ! write(message,'(a,i0,a,i0,a)') trim(message)//'mismatch for HRU ', idStruct%var(iLookID%hruId), '(param HRU = ', hruId(indxHRU), ')' + ! err=20; return + ! endif + + ! error check + ! else + ! err = 20; message = 'run mode not recognized'; return; + ! end if + + endif ! if the HRU id + + end do ! looping through variables in the file - endif ! reading the basin parameters + ! ********************************************************************************************** + ! * read the local parameters and the basin parameters + ! ********************************************************************************************** - end do ! (looping through the parameters in the NetCDF file) + ! loop through the parameters in the NetCDF file + do ivarid=1,nVars + + ! get the parameter name + err=nf90_inquire_variable(ncid, ivarid, name=parName) + call netcdf_err(err,message); if (err/=0) then; err=20; return; end if + + ! get the local parameters + ixParam = get_ixparam( trim(parName) ) + if(ixParam/=integerMissing)then + + ! ********************************************************************************************** + ! * read the local parameters + ! ********************************************************************************************** + + ! get the variable shape + err=nf90_inquire_variable(ncid, ivarid, nDims=nDims, dimids=idim_list) + if(err/=0)then; message=trim(message)//trim(cmessage); return; end if + + ! get the length of the depth dimension (if it exists) + if(nDims==2)then + + ! get the information on the 2nd dimension for 2-d variables + err=nf90_inquire_dimension(ncid, idim_list(2), dimName, nSoil_file) + if(err/=0)then; message=trim(message)//trim(cmessage); return; end if + + ! check that it is the depth dimension + if(trim(dimName)/='depth')then + message=trim(message)//'expect 2nd dimension of 2-d variable to be depth (dimension name = '//trim(dimName)//')' + err=20; return + endif + + ! check that the dimension length is correct + if(size(mparStruct%var(ixParam)%dat) /= nSoil_file)then + message=trim(message)//'unexpected number of soil layers in parameter file' + err=20; return + endif + + ! define parameter length + parLength = nSoil_file + + else + parLength = 1 + endif ! if two dimensions + + ! allocate space for model parameters + allocate(parVector(parLength),stat=err) + if(err/=0)then + message=trim(message)//'problem allocating space for parameter vector' + err=20; return + endif + + ! map to the GRUs and HRUs + !iGRU=index_map(indxHRU)%gru_ix + localHRU_ix=index_map(indxHRU)%localHRU_ix + fHRU = gru_struc(indxGRU)%hruInfo(localHRU_ix)%hru_nc + + ! read parameter data + select case(nDims) + case(1); err=nf90_get_var(ncid, ivarid, parVector, start=(/fHRU/), count=(/1/) ) + case(2); err=nf90_get_var(ncid, ivarid, parVector, start=(/fHRU,1/), count=(/1,nSoil_file/) ) + case default; err=20; message=trim(message)//'unexpected number of dimensions for parameter '//trim(parName) + end select + + ! error check for the parameter read + if(err/=0)then; message=trim(message)//trim(cmessage); return; end if + + ! populate parameter structures + select case(nDims) + case(1); mparStruct%var(ixParam)%dat(:) = parVector(1) ! also distributes scalar across depth dimension + case(2); mparStruct%var(ixParam)%dat(:) = parVector(:) + case default; err=20; message=trim(message)//'unexpected number of dimensions for parameter '//trim(parName) + end select + + ! end do ! looping through HRUs + + ! deallocate space for model parameters + deallocate(parVector,stat=err) + if(err/=0)then + message=trim(message)//'problem deallocating space for parameter vector' + err=20; return + endif - ! Now we must close the netcdf file - call nc_file_close(ncid,err,message) - if(err/=0)then;message=trim(message)//trim(cmessage);return;end if + ! ********************************************************************************************** + ! * read the basin parameters + ! ********************************************************************************************** + + ! get the basin parameters + else - end subroutine read_param4chm + ! get the parameter index + ixParam = get_ixbpar( trim(parName) ) + + ! allow extra variables in the file that are not used + if(ixParam==integerMissing) cycle + + ! allocate space for model parameters + allocate(parVector(nGRU_file),stat=err) + if(err/=0)then + message=trim(message)//'problem allocating space for parameter vector' + err=20; return + endif + + ! read parameter data + err=nf90_get_var(ncid, ivarid, parVector ) + if(err/=0)then; message=trim(message)//trim(cmessage); return; end if + + ! populate parameter structures + if (iRunMode==iRunModeGRU) then + !do iGRU=1,nGRU + bparStruct%var(ixParam) = parVector(indxGRU+startGRU-1) + !end do ! looping through GRUs + else if (iRunMode==iRunModeFull) then + !do iGRU=1,nGRU + bparStruct%var(ixParam) = parVector(indxGRU) + !end do ! looping through GRUs + else if (iRunMode==iRunModeHRU) then + err = 20; message='checkHRU run mode not working'; return; + endif + + ! deallocate space for model parameters + deallocate(parVector,stat=err) + if(err/=0)then + message=trim(message)//'problem deallocating space for parameter vector' + err=20; return + endif + + endif ! reading the basin parameters + + end do ! (looping through the parameters in the NetCDF file) + + ! Now we must close the netcdf file + call nc_file_close(ncid,err,message) + if(err/=0)then;message=trim(message)//trim(cmessage);return;end if + + end subroutine read_param end module read_param4chm_module -- GitLab