From ab093b892c6c97faea65a74335a500b034801152 Mon Sep 17 00:00:00 2001 From: Andy Wood Date: Mon, 21 Sep 2020 20:57:11 -0600 Subject: [PATCH 01/13] Fix restart file write/read to include basin hillslope runoff variables, allowing for initialized predictions that use basin runoff. Prior statefiles will no longer match this version, but also won't support exact restarts for basin variables. --- build/source/driver/summa_alarms.f90 | 2 +- build/source/driver/summa_init.f90 | 2 +- build/source/driver/summa_restart.f90 | 11 +-- build/source/driver/summa_writeOutput.f90 | 3 +- build/source/dshare/get_ixname.f90 | 2 +- build/source/engine/allocspace.f90 | 6 +- build/source/engine/updateVars.f90 | 16 ++++ build/source/netcdf/modelwrite.f90 | 67 +++++++++++----- build/source/netcdf/read_icond.f90 | 97 ++++++++++++++++++++--- 9 files changed, 166 insertions(+), 40 deletions(-) diff --git a/build/source/driver/summa_alarms.f90 b/build/source/driver/summa_alarms.f90 index 0f70ffed8..0d9c27864 100755 --- a/build/source/driver/summa_alarms.f90 +++ b/build/source/driver/summa_alarms.f90 @@ -123,7 +123,7 @@ subroutine summa_setWriteAlarms(oldTime, newTime, endTime, & ! time vect case(ixRestart_end); printRestart = (newTime(iLookTIME%im) == endTime(iLookTIME%im) .and. & newTime(iLookTIME%id) == endTime(iLookTIME%id) .and. & newTime(iLookTIME%ih) == endTime(iLookTIME%ih) .and. & - newTime(iLookTIME%imin) == endTime(iLookTIME%imin)) + newTime(iLookTIME%imin) == endTime(iLookTIME%imin)) ! newTime does not have a '24h', won't write ending state if end_h=24 case(ixRestart_never); printRestart = .false. case default; err=20; message=trim(message)//'unable to identify option for the restart file'; return end select diff --git a/build/source/driver/summa_init.f90 b/build/source/driver/summa_init.f90 index 7232288fc..1f65af736 100755 --- a/build/source/driver/summa_init.f90 +++ b/build/source/driver/summa_init.f90 @@ -121,7 +121,7 @@ subroutine summa_initialize(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 variables ! primary data structures (scalars) timeStruct => summa1_struc%timeStruct , & ! x%var(:) -- model time data diff --git a/build/source/driver/summa_restart.f90 b/build/source/driver/summa_restart.f90 index f179603a0..5080ec33c 100755 --- a/build/source/driver/summa_restart.f90 +++ b/build/source/driver/summa_restart.f90 @@ -107,7 +107,7 @@ subroutine summa_readRestart(summa1_struc, err, message) ! identify the start of the writing call date_and_time(values=startRestart) - + ! ***************************************************************************** ! *** read/check initial conditions ! ***************************************************************************** @@ -116,7 +116,7 @@ subroutine summa_readRestart(summa1_struc, err, message) if(STATE_PATH == '') then restartFile = trim(SETTINGS_PATH)//trim(MODEL_INITCOND) else - restartFile = trim(STATE_PATH)//trim(MODEL_INITCOND) + restartFile = trim(STATE_PATH)//trim(MODEL_INITCOND) endif ! read initial conditions @@ -124,6 +124,7 @@ subroutine summa_readRestart(summa1_struc, err, message) nGRU, & ! intent(in): number of response units mparStruct, & ! intent(in): model parameters progStruct, & ! intent(inout): model prognostic variables + bvarStruct, & ! intent(inout): model basin (GRU) variables indxStruct, & ! intent(inout): model indices err,cmessage) ! intent(out): error control if(err/=0)then; message=trim(message)//trim(cmessage); return; endif @@ -143,7 +144,7 @@ subroutine summa_readRestart(summa1_struc, err, message) ! *** compute ancillary variables ! ***************************************************************************** - ! loop through local HRUs + ! loop through HRUs do iHRU=1,gru_struc(iGRU)%hruCount ! re-calculate height of each layer @@ -178,7 +179,7 @@ subroutine summa_readRestart(summa1_struc, err, message) ! NOTE: canopy drip from the previous time step is used to compute throughfall for the current time step fluxStruct%gru(iGRU)%hru(iHRU)%var(iLookFLUX%scalarCanopyLiqDrainage)%dat(1) = 0._dp ! not used - end do ! (looping through HRUs) + end do ! end looping through HRUs ! ***************************************************************************** ! *** initialize aquifer storage @@ -225,7 +226,7 @@ subroutine summa_readRestart(summa1_struc, err, message) dt_init%gru(iGRU)%hru(iHRU) = progStruct%gru(iGRU)%hru(iHRU)%var(iLookPROG%dt_init)%dat(1) ! seconds end do - end do ! (looping through GRUs) + end do ! end looping through GRUs ! ***************************************************************************** ! *** finalize diff --git a/build/source/driver/summa_writeOutput.f90 b/build/source/driver/summa_writeOutput.f90 index 8dbeff23e..121fde7bc 100755 --- a/build/source/driver/summa_writeOutput.f90 +++ b/build/source/driver/summa_writeOutput.f90 @@ -287,7 +287,8 @@ subroutine summa_writeOutputFiles(modelTimeStep, summa1_struc, err, message) restartFile=trim(STATE_PATH)//trim(OUTPUT_PREFIX)//'_restart_'//trim(timeString)//trim(output_fileSuffix)//'.nc' endif - call writeRestart(restartFile,nGRU,nHRU,prog_meta,progStruct,maxLayers,maxSnowLayers,indx_meta,indxStruct,err,cmessage) + !call writeRestart(restartFile,nGRU,nHRU,prog_meta,progStruct,maxLayers,maxSnowLayers,indx_meta,indxStruct,err,cmessage) + call writeRestart(restartFile,nGRU,nHRU,prog_meta,progStruct,bvar_meta,bvarStruct,maxLayers,maxSnowLayers,indx_meta,indxStruct,err,cmessage) if(err/=0)then; message=trim(message)//trim(cmessage); return; endif end if diff --git a/build/source/dshare/get_ixname.f90 b/build/source/dshare/get_ixname.f90 index 926404ad9..1d136c8ac 100755 --- a/build/source/dshare/get_ixname.f90 +++ b/build/source/dshare/get_ixname.f90 @@ -977,7 +977,7 @@ subroutine get_ixUnknown(varName,typeName,vDex,err,message) message='get_ixUnknown/' ! loop through all structure types to find the one with the given variable name - ! pill variable index plus return which structure it was found in + ! poll variable index plus return which structure it was found in do iStruc = 1,size(structInfo) select case(trim(structInfo(iStruc)%structName)) case ('time' ); vDex = get_ixTime(trim(varName)) diff --git a/build/source/engine/allocspace.f90 b/build/source/engine/allocspace.f90 index e7e41ed5d..96ae8b5fd 100755 --- a/build/source/engine/allocspace.f90 +++ b/build/source/engine/allocspace.f90 @@ -66,9 +66,9 @@ module allocspace_module public::allocLocal public::resizeData -! define fixed dimensions -integer(i4b),parameter :: nBand=2 ! number of spectral bands -integer(i4b),parameter :: nTimeDelay=2000 ! number of elements in the time delay histogram +! define fixed dimensions (move these to globalData module?) +integer(i4b),parameter,public :: nBand=2 ! number of spectral bands +integer(i4b),parameter,public :: nTimeDelay=2000 ! number of elements in the time delay histogram ! ----------------------------------------------------------------------------------------------------------------------------------- contains diff --git a/build/source/engine/updateVars.f90 b/build/source/engine/updateVars.f90 index adee10fd6..c024f1cd2 100755 --- a/build/source/engine/updateVars.f90 +++ b/build/source/engine/updateVars.f90 @@ -713,6 +713,22 @@ subroutine xTempSolve(& ! subroutine starts here residual = -heatCap*(xTemp - tempInit) + meltNrg*(volFracIceTrial - volFracIceInit) ! J m-3 derivative = heatCap + LH_fus*iden_water*dLiq_dT ! J m-3 K-1 + + ! check validity of residual ... + ! informational only: if nan, the sim will start to error out from calling routine + if( ieee_is_nan(residual) )then + print*, '--------' + print*, 'ERROR: residual is not valid in xTempSolve' + print*, 'heatCap', heatCap + print*, 'xTemp', xTemp + print*, 'tempInit', tempInit + print*, 'meltNrg', meltNrg + print*, 'volFracIceTrial', volFracIceTrial + print*, 'volFracIceInit', volFracIceInit + print*, 'dLiq_dT', dLiq_dT + print*, '--------' + endif + end subroutine xTempSolve end module updateVars_module diff --git a/build/source/netcdf/modelwrite.f90 b/build/source/netcdf/modelwrite.f90 index dc417a9f4..665c39c91 100755 --- a/build/source/netcdf/modelwrite.f90 +++ b/build/source/netcdf/modelwrite.f90 @@ -356,7 +356,7 @@ subroutine writeBasin(iGRU,finalizeStats,outputTimestep,meta,stat,dat,map,err,me ! check that the variable is desired if (iStat==integerMissing.or.trim(meta(iVar)%varName)=='unknown') cycle - ! stats/dats output - select data type + ! stats/data output - select data type select case (meta(iVar)%varType) case (iLookVarType%scalarv) @@ -440,6 +440,8 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file nHRU, & ! intent(in): number of HRUs prog_meta, & ! intent(in): prognostics metadata prog_data, & ! intent(in): prognostics data + bvar_meta, & ! intent(in): basin (gru) variable metadata + bvar_data, & ! intent(in): basin (gru) variable data maxLayers, & ! intent(in): maximum number of layers maxSnowLayers, & ! intent(in): maximum number of snow layers indx_meta, & ! intent(in): index metadata @@ -452,19 +454,23 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file ! access named variables defining elements in the data structures USE var_lookup,only:iLookINDEX ! named variables for structure elements USE var_lookup,only:iLookVarType ! named variables for structure elements + USE var_lookup,only:iLookBVAR ! named variables for structure elements ! constants USE globalData,only:gru_struc ! gru-hru mapping structures ! external routines USE netcdf_util_module,only:nc_file_close ! close netcdf file USE netcdf_util_module,only:nc_file_open ! open netcdf file + !USE allocGlobal,only:nTimeDelay ! parameter setting for length of overland routing histogram implicit none ! -------------------------------------------------------------------------------------------------------- ! input character(len=256),intent(in) :: filename ! name of the restart file integer(i4b),intent(in) :: nGRU ! number of GRUs integer(i4b),intent(in) :: nHRU ! number of HRUs - type(var_info),intent(in) :: prog_meta(:) ! metadata + type(var_info),intent(in) :: prog_meta(:) ! prognostic variable metadata type(gru_hru_doubleVec),intent(in) :: prog_data ! prognostic vars + type(var_info),intent(in) :: bvar_meta(:) ! basin variable metadata + type(gru_doubleVec),intent(in) :: bvar_data ! basin variables type(var_info),intent(in) :: indx_meta(:) ! metadata type(gru_hru_intVec),intent(in) :: indx_data ! indexing vars ! output: error control @@ -486,10 +492,14 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file integer(i4b) :: maxSnow ! maximum number of snow layers integer(i4b) :: maxSoil ! maximum number of soil layers integer(i4b) :: nLayers ! number of total layers + integer(i4b) :: nTDH ! number of points in time-delay histogram integer(i4b),parameter :: nSpectral=2 ! number of spectal bands integer(i4b),parameter :: nScalar=1 ! size of a scalar + integer(i4b) :: nProgVars ! number of prognostic variables written to state file integer(i4b) :: hruDimID ! variable dimension ID + integer(i4b) :: gruDimID ! variable dimension ID + integer(i4b) :: tdhDimID ! variable dimension ID integer(i4b) :: scalDimID ! variable dimension ID integer(i4b) :: specDimID ! variable dimension ID integer(i4b) :: midSnowDimID ! variable dimension ID @@ -500,6 +510,8 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file integer(i4b) :: ifcTotoDimID ! variable dimension ID character(len=32),parameter :: hruDimName ='hru' ! dimension name for HRUs + character(len=32),parameter :: gruDimName ='gru' ! dimension name for GRUs + character(len=32),parameter :: tdhDimName ='tdh' ! dimension name for time-delay basin variables character(len=32),parameter :: scalDimName ='scalarv' ! dimension name for scalar data character(len=32),parameter :: specDimName ='spectral' ! dimension name for spectral bands character(len=32),parameter :: midSnowDimName='midSnow' ! dimension name for snow-only layers @@ -507,7 +519,7 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file character(len=32),parameter :: midTotoDimName='midToto' ! dimension name for layered varaiables character(len=32),parameter :: ifcSnowDimName='ifcSnow' ! dimension name for snow-only layers character(len=32),parameter :: ifcSoilDimName='ifcSoil' ! dimension name for soil-only layers - character(len=32),parameter :: ifcTotoDimName='ifcToto' ! dimension name for layered varaiables + character(len=32),parameter :: ifcTotoDimName='ifcToto' ! dimension name for layered variables integer(i4b) :: cHRU ! count of HRUs integer(i4b) :: iHRU ! index of HRUs @@ -520,34 +532,40 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file ! initialize error control err=0; message='writeRestart/' - ! size of prog vector - allocate(ncVarID(size(prog_meta))) + ! size of prognostic variable vector + nProgVars = size(prog_meta) + allocate(ncVarID(nProgVars+2)) ! include 2 additional basin variables in ID array ! maximum number of soil layers maxSoil = gru_struc(1)%hruInfo(1)%nSoil ! maximum number of snow layers maxSnow = maxSnowLayers + + ! length of time delay histogram + nTDH = 2000 ! figure a way to get this param from allocGlobal into here ! create file err = nf90_create(trim(filename),nf90_classic_model,ncid) message='iCreate[create]'; call netcdf_err(err,message); if(err/=0)return ! define dimensions - err = nf90_def_dim(ncid,trim(hruDimName) ,nHRU , hruDimID) ; message='iCreate[hru]' ;call netcdf_err(err,message); if(err/=0)return - err = nf90_def_dim(ncid,trim(scalDimName) ,nScalar , scalDimID); message='iCreate[scalar]' ;call netcdf_err(err,message); if(err/=0)return - err = nf90_def_dim(ncid,trim(specDimName) ,nSpectral , specDimID); message='iCreate[spectral]';call netcdf_err(err,message); if(err/=0)return - err = nf90_def_dim(ncid,trim(midSoilDimName),maxSoil ,midSoilDimID); message='iCreate[ifcSoil]' ;call netcdf_err(err,message); if(err/=0)return - err = nf90_def_dim(ncid,trim(midTotoDimName),maxLayers ,midTotoDimID); message='iCreate[midToto]' ;call netcdf_err(err,message); if(err/=0)return - err = nf90_def_dim(ncid,trim(ifcSoilDimName),maxSoil+1 ,ifcSoilDimID); message='iCreate[ifcSoil]' ;call netcdf_err(err,message); if(err/=0)return - err = nf90_def_dim(ncid,trim(ifcTotoDimName),maxLayers+1,ifcTotoDimID); message='iCreate[ifcToto]' ;call netcdf_err(err,message); if(err/=0)return - if (maxSnow>0) err = nf90_def_dim(ncid,trim(midSnowDimName),maxSnow ,midSnowDimID); message='iCreate[ifcSnow]' ;call netcdf_err(err,message); if(err/=0)return - if (maxSnow>0) err = nf90_def_dim(ncid,trim(ifcSnowDimName),maxSnow+1 ,ifcSnowDimID); message='iCreate[ifcSnow]' ;call netcdf_err(err,message); if(err/=0)return + err = nf90_def_dim(ncid,trim(hruDimName) ,nHRU , hruDimID); message='iCreate[hru]' ; call netcdf_err(err,message); if(err/=0)return + err = nf90_def_dim(ncid,trim(gruDimName) ,nGRU , gruDimID); message='iCreate[gru]' ; call netcdf_err(err,message); if(err/=0)return + err = nf90_def_dim(ncid,trim(tdhDimName) ,nTDH , tdhDimID); message='iCreate[tdh]' ; call netcdf_err(err,message); if(err/=0)return + err = nf90_def_dim(ncid,trim(scalDimName) ,nScalar , scalDimID); message='iCreate[scalar]' ; call netcdf_err(err,message); if(err/=0)return + err = nf90_def_dim(ncid,trim(specDimName) ,nSpectral , specDimID); message='iCreate[spectral]'; call netcdf_err(err,message); if(err/=0)return + err = nf90_def_dim(ncid,trim(midSoilDimName),maxSoil ,midSoilDimID); message='iCreate[ifcSoil]' ; call netcdf_err(err,message); if(err/=0)return + err = nf90_def_dim(ncid,trim(midTotoDimName),maxLayers ,midTotoDimID); message='iCreate[midToto]' ; call netcdf_err(err,message); if(err/=0)return + err = nf90_def_dim(ncid,trim(ifcSoilDimName),maxSoil+1 ,ifcSoilDimID); message='iCreate[ifcSoil]' ; call netcdf_err(err,message); if(err/=0)return + err = nf90_def_dim(ncid,trim(ifcTotoDimName),maxLayers+1,ifcTotoDimID); message='iCreate[ifcToto]' ; call netcdf_err(err,message); if(err/=0)return + if (maxSnow>0) err = nf90_def_dim(ncid,trim(midSnowDimName),maxSnow ,midSnowDimID); message='iCreate[ifcSnow]' ; call netcdf_err(err,message); if(err/=0)return + if (maxSnow>0) err = nf90_def_dim(ncid,trim(ifcSnowDimName),maxSnow+1 ,ifcSnowDimID); message='iCreate[ifcSnow]' ; call netcdf_err(err,message); if(err/=0)return ! re-initialize error control err=0; message='writeRestart/' ! define prognostic variables - do iVar = 1,size(prog_meta) + do iVar = 1,nProgVars if (prog_meta(iVar)%varType==iLookvarType%unknown) cycle ! define variable @@ -577,6 +595,14 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file call netcdf_err(err,message) end do ! iVar + + ! define selected basin variables (derived) -- e.g., hillslope routing + err = nf90_def_var(ncid, trim(bvar_meta(iLookBVAR%routingRunoffFuture)%varName), nf90_double, (/gruDimID, tdhDimID /), ncVarID(nProgVars+1)) + err = nf90_put_att(ncid,ncVarID(nProgVars+1),'long_name',trim(bvar_meta(iLookBVAR%routingRunoffFuture)%vardesc)); call netcdf_err(err,message) + err = nf90_put_att(ncid,ncVarID(nProgVars+1),'units' ,trim(bvar_meta(iLookBVAR%routingRunoffFuture)%varunit)); call netcdf_err(err,message) + err = nf90_def_var(ncid, trim(bvar_meta(iLookBVAR%routingFractionFuture)%varName), nf90_double, (/gruDimID, tdhDimID /), ncVarID(nProgVars+2)) + err = nf90_put_att(ncid,ncVarID(nProgVars+2),'long_name',trim(bvar_meta(iLookBVAR%routingFractionFuture)%vardesc)); call netcdf_err(err,message) + err = nf90_put_att(ncid,ncVarID(nProgVars+2),'units' ,trim(bvar_meta(iLookBVAR%routingFractionFuture)%varunit)); call netcdf_err(err,message) ! define index variables - snow err = nf90_def_var(ncid,trim(indx_meta(iLookIndex%nSnow)%varName),nf90_int,(/hruDimID/),ncSnowID); call netcdf_err(err,message) @@ -644,14 +670,19 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file call netcdf_err(err,message); if (err/=0) return err=0; message='writeRestart/' - end do ! iVar + end do ! iVar loop ! write index variables err=nf90_put_var(ncid,ncSnowID,(/indx_data%gru(iGRU)%hru(iHRU)%var(iLookIndex%nSnow)%dat/),start=(/cHRU/),count=(/1/)) err=nf90_put_var(ncid,ncSoilID,(/indx_data%gru(iGRU)%hru(iHRU)%var(iLookIndex%nSoil)%dat/),start=(/cHRU/),count=(/1/)) - end do ! iGRU - end do ! iHRU + end do ! iHRU loop + + ! write selected basin variables + err=nf90_put_var(ncid,ncVarID(nProgVars+1),(/bvar_data%gru(iGRU)%var(iLookBVAR%routingRunoffFuture)%dat/), start=(/iGRU/),count=(/1,nTDH/)) + err=nf90_put_var(ncid,ncVarID(nProgVars+2),(/bvar_data%gru(iGRU)%var(iLookBVAR%routingFractionFuture)%dat/),start=(/iGRU/),count=(/1,nTDH/)) + + end do ! iGRU loop ! close file call nc_file_close(ncid,err,cmessage) diff --git a/build/source/netcdf/read_icond.f90 b/build/source/netcdf/read_icond.f90 index 6559af4a7..39267112a 100755 --- a/build/source/netcdf/read_icond.f90 +++ b/build/source/netcdf/read_icond.f90 @@ -82,7 +82,7 @@ subroutine read_icond_nlayers(iconFile,nGRU,indx_meta,err,message) err = nf90_inq_dimid(ncID,"hru",dimId); if(err/=nf90_noerr)then; message=trim(message)//'problem finding hru dimension/'//trim(nf90_strerror(err)); return; end if err = nf90_inquire_dimension(ncID,dimId,len=fileHRU); if(err/=nf90_noerr)then; message=trim(message)//'problem reading hru dimension/'//trim(nf90_strerror(err)); return; end if - ! allocate sotrage for reading from file + ! allocate storage for reading from file allocate(snowData(fileHRU)) allocate(soilData(fileHRU)) snowData = 0 @@ -144,6 +144,7 @@ subroutine read_icond(iconFile, & ! intent(in): name of nGRU, & ! intent(in): number of GRUs mparData, & ! intent(in): model parameters progData, & ! intent(inout): model prognostic variables + bvarData, & ! intent(inout): model basin (GRU) variables indxData, & ! intent(inout): model indices err,message) ! intent(out): error control ! -------------------------------------------------------------------------------------------------------- @@ -152,8 +153,10 @@ subroutine read_icond(iconFile, & ! intent(in): name of USE var_lookup,only:iLookVarType ! variable lookup structure USE var_lookup,only:iLookPROG ! variable lookup structure USE var_lookup,only:iLookPARAM ! variable lookup structure + USE var_lookup,only:iLookBVAR ! variable lookup structure USE var_lookup,only:iLookINDEX ! variable lookup structure USE globalData,only:prog_meta ! metadata for prognostic variables + USE globalData,only:bvar_meta ! metadata for basin (GRU) variables USE globalData,only:gru_struc ! gru-hru mapping structures USE globaldata,only:iname_soil,iname_snow ! named variables to describe the type of layer USE netcdf_util_module,only:nc_file_open ! open netcdf file @@ -161,6 +164,7 @@ subroutine read_icond(iconFile, & ! intent(in): name of USE netcdf_util_module,only:netcdf_err ! netcdf error handling USE data_types,only:gru_hru_doubleVec ! full double precision structure USE data_types,only:gru_hru_intVec ! full integer structure + USE data_types,only:gru_doubleVec ! gru-length double precision structure (basin variables) USE data_types,only:var_dlength ! double precision structure for a single HRU USE data_types,only:var_info ! metadata USE get_ixName_module,only:get_varTypeName ! to access type strings for error messages @@ -174,6 +178,7 @@ subroutine read_icond(iconFile, & ! intent(in): name of integer(i4b) ,intent(in) :: nGRU ! number of grouped response units in simulation domain type(gru_hru_doubleVec),intent(in) :: mparData ! model parameters type(gru_hru_doubleVec),intent(inout) :: progData ! model prognostic variables + type(gru_doubleVec) ,intent(inout) :: bvarData ! model basin (GRU) variables type(gru_hru_intVec) ,intent(inout) :: indxData ! model indices integer(i4b) ,intent(out) :: err ! error code character(*) ,intent(out) :: message ! returned error message @@ -181,9 +186,12 @@ subroutine read_icond(iconFile, & ! intent(in): name of ! locals character(len=256) :: cmessage ! downstream error message integer(i4b) :: fileHRU ! number of HRUs in file - integer(i4b) :: iVar ! loop index + integer(i4b) :: fileGRU ! number of GRUs in file + integer(i4b) :: iVar, i ! loop indices + integer(i4b),dimension(2) :: ndx ! intermediate array of loop indices integer(i4b) :: iGRU ! loop index integer(i4b) :: iHRU ! loop index + integer(i4b) :: iTDH ! loop index (maybe not needed) integer(i4b) :: dimID ! varible dimension ids integer(i4b) :: ncVarID ! variable ID in netcdf file character(256) :: dimName ! not used except as a placeholder in call to inq_dim function @@ -194,18 +202,25 @@ subroutine read_icond(iconFile, & ! intent(in): name of integer(i4b) :: iHRU_global ! index of HRU in the netcdf file real(dp),allocatable :: varData(:,:) ! variable data storage integer(i4b) :: nSoil, nSnow, nToto ! # layers + integer(i4b) :: nTDH ! number of points in time-delay histogram integer(i4b) :: iLayer,jLayer ! layer indices integer(i4b),parameter :: nBand=2 ! number of spectral bands + integer(i4b) :: nProgVars ! number of prognostic variables written to state file character(len=32),parameter :: scalDimName ='scalarv' ! dimension name for scalar data character(len=32),parameter :: midSoilDimName='midSoil' ! dimension name for soil-only layers character(len=32),parameter :: midTotoDimName='midToto' ! dimension name for layered varaiables character(len=32),parameter :: ifcTotoDimName='ifcToto' ! dimension name for layered varaiables + character(len=32),parameter :: tdhDimName ='tdh' ! dimension name for time-delay basin variables (maybe not needed) + ! -------------------------------------------------------------------------------------------------------- ! Start procedure here err=0; message="read_icond/" + + ! set length of time delay histogram + nTDH = 2000 ! TODO figure a way to get this param from allocGlobal into here ! -------------------------------------------------------------------------------------------------------- ! (1) read the file @@ -218,6 +233,15 @@ subroutine read_icond(iconFile, & ! intent(in): name of err = nf90_inq_dimid(ncID,"hru",dimID); if(err/=nf90_noerr)then; message=trim(message)//'problem finding hru dimension/'//trim(nf90_strerror(err)); return; end if err = nf90_inquire_dimension(ncID,dimID,len=fileHRU); if(err/=nf90_noerr)then; message=trim(message)//'problem reading hru dimension/'//trim(nf90_strerror(err)); return; end if + ! get number of GRUs in file + err = nf90_inq_dimid(ncID,"gru",dimID); if(err/=nf90_noerr)then; message=trim(message)//'problem finding gru dimension/'//trim(nf90_strerror(err)); return; end if + err = nf90_inquire_dimension(ncID,dimID,len=fileGRU); if(err/=nf90_noerr)then; message=trim(message)//'problem reading gru dimension/'//trim(nf90_strerror(err)); return; end if + + ! get dimension of time delay histogram (TDH) in file + err = nf90_inq_dimid(ncID,"tdh",dimID); if(err/=nf90_noerr)then; message=trim(message)//'problem finding tdh dimension/'//trim(nf90_strerror(err)); return; end if + err = nf90_inquire_dimension(ncID,dimID,len=nTDH); if(err/=nf90_noerr)then; message=trim(message)//'problem reading tdh dimension/'//trim(nf90_strerror(err)); return; end if + ! need check via hardwired value elsewhere in code + ! loop through prognostic variables do iVar = 1,size(prog_meta) @@ -256,9 +280,9 @@ subroutine read_icond(iconFile, & ! intent(in): name of err = nf90_inquire_dimension(ncID,dimID,dimName,dimLen); call netcdf_err(err,message) if(err/=0)then; message=trim(message)//': problem getting the dimension length'; return; endif - ! iniitialize the variable data + ! initialize the variable data allocate(varData(fileHRU,dimLen),stat=err) - if(err/=0)then; message=trim(message)//'problem allocating variable data'; return; endif + if(err/=0)then; message=trim(message)//'problem allocating HRU variable data'; return; endif ! get data err = nf90_get_var(ncID,ncVarID,varData); call netcdf_err(err,message) @@ -269,8 +293,9 @@ subroutine read_icond(iconFile, & ! intent(in): name of do iGRU = 1,nGRU do iHRU = 1,gru_struc(iGRU)%hruCount - iHRU_global = gru_struc(iGRU)%hruInfo(iHRU)%hru_nc - iHRU_local = (iHRU_global - ixHRUfile_min) + 1 + iHRU_global = gru_struc(iGRU)%hruInfo(iHRU)%hru_nc + iHRU_local = (iHRU_global - ixHRUfile_min) + 1 + ! get the number of layers nSnow = gru_struc(iGRU)%hruInfo(iHRU)%nSnow nSoil = gru_struc(iGRU)%hruInfo(iHRU)%nSoil @@ -319,10 +344,10 @@ subroutine read_icond(iconFile, & ! intent(in): name of ! deallocate storage vector for next variable deallocate(varData, stat=err) - if(err/=0)then; message=trim(message)//'problem deallocating variable data'; return; endif - - end do ! iVar + if(err/=0)then; message=trim(message)//'problem deallocating HRU variable data'; return; endif + end do ! end looping through prognostic variables (iVar) + ! -------------------------------------------------------------------------------------------------------- ! (2) set number of layers ! -------------------------------------------------------------------------------------------------------- @@ -342,7 +367,7 @@ subroutine read_icond(iconFile, & ! intent(in): name of end do ! -------------------------------------------------------------------------------------------------------- - ! (3) update soil layers + ! (3) update soil layers (diagnostic variables) ! -------------------------------------------------------------------------------------------------------- ! loop through GRUs and HRUs @@ -376,6 +401,58 @@ subroutine read_icond(iconFile, & ! intent(in): name of end do ! looping throuygh HRUs end do ! looping through GRUs + + ! -------------------------------------------------------------------------------------------------------- + ! (2) now get the two basin runoff variables + ! -------------------------------------------------------------------------------------------------------- + + ! get the index in the file: single HRU + if(restartFileType/=singleHRU)then + + ! loop through specific basin variables + ndx = (/iLookBVAR%routingRunoffFuture, iLookBVAR%routingFractionFuture/) ! array of desired variable indices + do i = 1,size(ndx) + iVar = ndx(i) + + ! get tdh dimension Id in file (should be 'tdh') + err = nf90_inq_dimid(ncID,trim(tdhDimName), dimID); + if(err/=0)then; message=trim(message)//': problem with dimension ids for tdh vars'; return; endif + + ! get the tdh dimension length (dimName and dimLen are outputs of this call) + err = nf90_inquire_dimension(ncID,dimID,dimName,dimLen); call netcdf_err(err,message) + if(err/=0)then; message=trim(message)//': problem getting the dimension length for tdh vars'; return; endif + + ! get tdh-based variable id + err = nf90_inq_varid(ncID,trim(bvar_meta(iVar)%varName),ncVarID); call netcdf_err(err,message) + if(err/=0)then; message=trim(message)//': problem with getting basin variable id, var='//trim(bvar_meta(iVar)%varName); return; endif + + ! initialize the tdh variable data + allocate(varData(fileGRU,dimLen),stat=err) + if(err/=0)then; print*, 'err= ',err; message=trim(message)//'problem allocating GRU variable data'; return; endif + + ! get data + err = nf90_get_var(ncID,ncVarID,varData); call netcdf_err(err,message) + if(err/=0)then; message=trim(message)//': problem getting the data'; return; endif + + ! store data in basin var (bvar) structure + do iGRU = 1,nGRU + + ! put the data into data structures + bvarData%gru(iGRU)%var(iVar)%dat(1:nTDH) = varData(iGRU,1:nTDH) + ! check whether the first values is set to nf90_fill_double + if(any(abs(bvarData%gru(iGRU)%var(iVar)%dat(1:nTDH) - nf90_fill_double) < epsilon(varData)))then; err=20; endif + if(err==20)then; message=trim(message)//"data set to the fill value (name='"//trim(bvar_meta(iVar)%varName)//"')"; return; endif + + end do ! end iGRU loop + + ! deallocate temporary data array for next variable + deallocate(varData, stat=err) + if(err/=0)then; message=trim(message)//'problem deallocating GRU variable data'; return; endif + + end do ! end looping through basin variables + + endif ! end if case for not being a singleHRU run + end subroutine read_icond From 48ec3a32907f1ff099d7c1019831c213169d1c4c Mon Sep 17 00:00:00 2001 From: Andy Wood Date: Tue, 22 Sep 2020 08:35:17 -0600 Subject: [PATCH 02/13] increased tolerance on differences in height & layer thickness from 10e-12 which caused errors below float level precision in cold state files --- build/source/engine/check_icond.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/build/source/engine/check_icond.f90 b/build/source/engine/check_icond.f90 index ba61a20f1..7664270ff 100755 --- a/build/source/engine/check_icond.f90 +++ b/build/source/engine/check_icond.f90 @@ -272,8 +272,8 @@ subroutine check_icond(nGRU, & ! number of GRUs and HRU do iLayer=1,nLayers h1 = sum(progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%mLayerDepth)%dat(1:iLayer)) ! sum of the depths up to the current layer h2 = progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%iLayerHeight)%dat(iLayer) - progData%gru(iGRU)%hru(iHRU)%var(iLookPROG%iLayerHeight)%dat(0) ! difference between snow-atm interface and bottom of layer - if(abs(h1 - h2) > 1.e-12_dp)then - write(message,'(a,1x,i0)') trim(message)//'mis-match between layer depth and layer height [suggest round numbers in initial conditions file]; layer = ', iLayer + if(abs(h1 - h2) > 1.e-6_dp)then + write(message,'(a,1x,i0)') trim(message)//'mis-match between layer depth and layer height; layer = ', iLayer, '; sum depths = ',h1,'; height = ',h2 err=20; return end if end do From 6638488188210c3508ceeb5ed3c70b2be4eb3b08 Mon Sep 17 00:00:00 2001 From: Andy Wood Date: Tue, 22 Sep 2020 22:55:21 -0600 Subject: [PATCH 03/13] Refined state file basin variables fix: fraction future not needed in state file; moved nTimeDelay to globalData and set only once --- build/source/dshare/globalData.f90 | 4 ++++ build/source/engine/allocspace.f90 | 22 +++++++++++++--------- build/source/engine/check_icond.f90 | 1 + build/source/engine/var_derive.f90 | 10 +++++++--- build/source/netcdf/modelwrite.f90 | 19 ++++++------------- build/source/netcdf/read_icond.f90 | 22 ++++++++++++++-------- 6 files changed, 45 insertions(+), 33 deletions(-) diff --git a/build/source/dshare/globalData.f90 b/build/source/dshare/globalData.f90 index 5fb4301ac..7ff1a9ce7 100755 --- a/build/source/dshare/globalData.f90 +++ b/build/source/dshare/globalData.f90 @@ -336,5 +336,9 @@ MODULE globalData integer(i4b),parameter,public :: ncTime=1 ! time zone information from NetCDF file (timeOffset = longitude/15. - ncTimeOffset) integer(i4b),parameter,public :: utcTime=2 ! all times in UTC (timeOffset = longitude/15. hours) integer(i4b),parameter,public :: localTime=3 ! all times local (timeOffset = 0) + + ! define fixed dimensions (move these to globalData module?) + integer(i4b),parameter,public :: nBand=2 ! number of spectral bands + integer(i4b),parameter,public :: nTimeDelay=2000 ! number of hours in the time delay histogram (default: ~1 season = 24*365/4) END MODULE globalData diff --git a/build/source/engine/allocspace.f90 b/build/source/engine/allocspace.f90 index 96ae8b5fd..07a28d643 100755 --- a/build/source/engine/allocspace.f90 +++ b/build/source/engine/allocspace.f90 @@ -55,6 +55,9 @@ module allocspace_module USE globalData,only:integerMissing ! missing integer USE globalData,only:realMissing ! missing double precision number +USE globalData,only: nTimeDelay ! number of timesteps in the time delay histogram +USE globalData,only: nBand ! number of spectral bands + ! access variable types USE var_lookup,only:iLookVarType ! look up structure for variable typed USE var_lookup,only:maxvarFreq ! allocation dimension (output frequency) @@ -66,9 +69,6 @@ module allocspace_module public::allocLocal public::resizeData -! define fixed dimensions (move these to globalData module?) -integer(i4b),parameter,public :: nBand=2 ! number of spectral bands -integer(i4b),parameter,public :: nTimeDelay=2000 ! number of elements in the time delay histogram ! ----------------------------------------------------------------------------------------------------------------------------------- contains @@ -530,6 +530,7 @@ subroutine allocateDat_dp(metadata,nSnow,nSoil,nLayers, & ! input varData,err,message) ! output ! access subroutines USE get_ixName_module,only:get_varTypeName ! to access type strings for error messages + implicit none ! input variables type(var_info),intent(in) :: metadata(:) ! metadata structure @@ -543,12 +544,13 @@ subroutine allocateDat_dp(metadata,nSnow,nSoil,nLayers, & ! input ! local variables integer(i4b) :: iVar ! variable index integer(i4b) :: nVars ! number of variables in the metadata structure - ! initialize error control + +! initialize error control err=0; message='allocateDat_dp/' ! get the number of variables in the metadata structure nVars = size(metadata) - + ! loop through variables in the data structure do iVar=1,nVars @@ -608,13 +610,14 @@ subroutine allocateDat_int(metadata,nSnow,nSoil,nLayers, & ! input ! local variables integer(i4b) :: iVar ! variable index integer(i4b) :: nVars ! number of variables in the metadata structure - ! initialize error control + +! initialize error control err=0; message='allocateDat_int/' ! get the number of variables in the metadata structure nVars = size(metadata) - ! loop through variables in the data structure +! loop through variables in the data structure do iVar=1,nVars ! check allocated @@ -670,13 +673,14 @@ subroutine allocateDat_flag(metadata,nSnow,nSoil,nLayers, & ! input ! local variables integer(i4b) :: iVar ! variable index integer(i4b) :: nVars ! number of variables in the metadata structure - ! initialize error control + +! initialize error control err=0; message='allocateDat_flag/' ! get the number of variables in the metadata structure nVars = size(metadata) - ! loop through variables in the data structure +! loop through variables in the data structure do iVar=1,nVars ! check allocated diff --git a/build/source/engine/check_icond.f90 b/build/source/engine/check_icond.f90 index 7664270ff..9a1e9a779 100755 --- a/build/source/engine/check_icond.f90 +++ b/build/source/engine/check_icond.f90 @@ -207,6 +207,7 @@ subroutine check_icond(nGRU, & ! number of GRUs and HRU if(scalarTheta > theta_sat(iSoil)+xTol)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with total water fraction [liquid + ice] > theta_sat: layer = ',iLayer; err=20; return; end if case default + write(*,*) 'Cannot recognize case in initial vol water/ice check: type=', mlayerLayerType(iLayer) err=20; message=trim(message)//'cannot identify layer type'; return end select diff --git a/build/source/engine/var_derive.f90 b/build/source/engine/var_derive.f90 index cdf93f990..1b2fbcda2 100755 --- a/build/source/engine/var_derive.f90 +++ b/build/source/engine/var_derive.f90 @@ -381,6 +381,7 @@ end subroutine satHydCond subroutine fracFuture(bpar_data,bvar_data,err,message) ! external functions USE soil_utils_module,only:gammp ! compute the cumulative probabilty based on the Gamma distribution + implicit none ! input variables real(dp),intent(in) :: bpar_data(:) ! vector of basin-average model parameters @@ -397,7 +398,7 @@ subroutine fracFuture(bpar_data,bvar_data,err,message) real(dp) :: pSave ! cumulative probability at the start of the step real(dp) :: cumProb ! cumulative probability at the end of the step real(dp) :: sumFrac ! sum of runoff fractions in all steps - real(dp),parameter :: tolerFrac=0.01_dp ! tolerance for fractional runoff + real(dp),parameter :: tolerFrac=0.01_dp ! tolerance for missing fractional runoff by truncating histogram ! initialize error control err=0; message='fracFuture/' ! ---------------------------------------------------------------------------------- @@ -414,7 +415,7 @@ subroutine fracFuture(bpar_data,bvar_data,err,message) ! define time step dt = data_step ! length of the data step (s) - ! identify number of points in the time-delay histogram + ! identify number of points in the time-delay runoff variable (should be allocated match nTimeDelay) nTDH = size(runoffFuture) ! initialize runoffFuture @@ -451,10 +452,13 @@ subroutine fracFuture(bpar_data,bvar_data,err,message) exit end if end do ! (looping through future time steps) + ! check that we have enough bins sumFrac = sum(fractionFuture) if(abs(1._dp - sumFrac) > tolerFrac)then - message=trim(message)//'not enough bins for the time delay histogram -- fix hard-coded parameter in alloc_bvar' + write(*,*) 'fraction of basin runoff histogram being accounted for by time delay vector is ', sumFrac + write(*,*) 'this is less than allowed by tolerFrac = ', tolerFrac + message=trim(message)//'not enough bins for the time delay histogram -- fix hard-coded parameter in allocspace.f90' err=20; return end if ! ensure the fraction sums to one diff --git a/build/source/netcdf/modelwrite.f90 b/build/source/netcdf/modelwrite.f90 index 665c39c91..ee27a52af 100755 --- a/build/source/netcdf/modelwrite.f90 +++ b/build/source/netcdf/modelwrite.f90 @@ -28,7 +28,7 @@ module modelwrite_module USE nrtype ! missing values -USE globalData, only: integerMissing, realMissing +USE globalData,only: integerMissing, realMissing ! provide access to global data USE globalData,only:gru_struc ! gru->hru mapping structure @@ -460,7 +460,8 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file ! external routines USE netcdf_util_module,only:nc_file_close ! close netcdf file USE netcdf_util_module,only:nc_file_open ! open netcdf file - !USE allocGlobal,only:nTimeDelay ! parameter setting for length of overland routing histogram + USE globalData,only:nTimeDelay ! number of timesteps in the time delay histogram + implicit none ! -------------------------------------------------------------------------------------------------------- ! input @@ -492,7 +493,6 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file integer(i4b) :: maxSnow ! maximum number of snow layers integer(i4b) :: maxSoil ! maximum number of soil layers integer(i4b) :: nLayers ! number of total layers - integer(i4b) :: nTDH ! number of points in time-delay histogram integer(i4b),parameter :: nSpectral=2 ! number of spectal bands integer(i4b),parameter :: nScalar=1 ! size of a scalar integer(i4b) :: nProgVars ! number of prognostic variables written to state file @@ -534,7 +534,7 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file ! size of prognostic variable vector nProgVars = size(prog_meta) - allocate(ncVarID(nProgVars+2)) ! include 2 additional basin variables in ID array + allocate(ncVarID(nProgVars+1)) ! include 1 additional basin variable in ID array (possibly more later) ! maximum number of soil layers maxSoil = gru_struc(1)%hruInfo(1)%nSoil @@ -542,9 +542,6 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file ! maximum number of snow layers maxSnow = maxSnowLayers - ! length of time delay histogram - nTDH = 2000 ! figure a way to get this param from allocGlobal into here - ! create file err = nf90_create(trim(filename),nf90_classic_model,ncid) message='iCreate[create]'; call netcdf_err(err,message); if(err/=0)return @@ -552,7 +549,7 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file ! define dimensions err = nf90_def_dim(ncid,trim(hruDimName) ,nHRU , hruDimID); message='iCreate[hru]' ; call netcdf_err(err,message); if(err/=0)return err = nf90_def_dim(ncid,trim(gruDimName) ,nGRU , gruDimID); message='iCreate[gru]' ; call netcdf_err(err,message); if(err/=0)return - err = nf90_def_dim(ncid,trim(tdhDimName) ,nTDH , tdhDimID); message='iCreate[tdh]' ; call netcdf_err(err,message); if(err/=0)return + err = nf90_def_dim(ncid,trim(tdhDimName) ,nTimeDelay , tdhDimID); message='iCreate[tdh]' ; call netcdf_err(err,message); if(err/=0)return err = nf90_def_dim(ncid,trim(scalDimName) ,nScalar , scalDimID); message='iCreate[scalar]' ; call netcdf_err(err,message); if(err/=0)return err = nf90_def_dim(ncid,trim(specDimName) ,nSpectral , specDimID); message='iCreate[spectral]'; call netcdf_err(err,message); if(err/=0)return err = nf90_def_dim(ncid,trim(midSoilDimName),maxSoil ,midSoilDimID); message='iCreate[ifcSoil]' ; call netcdf_err(err,message); if(err/=0)return @@ -600,9 +597,6 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file err = nf90_def_var(ncid, trim(bvar_meta(iLookBVAR%routingRunoffFuture)%varName), nf90_double, (/gruDimID, tdhDimID /), ncVarID(nProgVars+1)) err = nf90_put_att(ncid,ncVarID(nProgVars+1),'long_name',trim(bvar_meta(iLookBVAR%routingRunoffFuture)%vardesc)); call netcdf_err(err,message) err = nf90_put_att(ncid,ncVarID(nProgVars+1),'units' ,trim(bvar_meta(iLookBVAR%routingRunoffFuture)%varunit)); call netcdf_err(err,message) - err = nf90_def_var(ncid, trim(bvar_meta(iLookBVAR%routingFractionFuture)%varName), nf90_double, (/gruDimID, tdhDimID /), ncVarID(nProgVars+2)) - err = nf90_put_att(ncid,ncVarID(nProgVars+2),'long_name',trim(bvar_meta(iLookBVAR%routingFractionFuture)%vardesc)); call netcdf_err(err,message) - err = nf90_put_att(ncid,ncVarID(nProgVars+2),'units' ,trim(bvar_meta(iLookBVAR%routingFractionFuture)%varunit)); call netcdf_err(err,message) ! define index variables - snow err = nf90_def_var(ncid,trim(indx_meta(iLookIndex%nSnow)%varName),nf90_int,(/hruDimID/),ncSnowID); call netcdf_err(err,message) @@ -679,8 +673,7 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file end do ! iHRU loop ! write selected basin variables - err=nf90_put_var(ncid,ncVarID(nProgVars+1),(/bvar_data%gru(iGRU)%var(iLookBVAR%routingRunoffFuture)%dat/), start=(/iGRU/),count=(/1,nTDH/)) - err=nf90_put_var(ncid,ncVarID(nProgVars+2),(/bvar_data%gru(iGRU)%var(iLookBVAR%routingFractionFuture)%dat/),start=(/iGRU/),count=(/1,nTDH/)) + err=nf90_put_var(ncid,ncVarID(nProgVars+1),(/bvar_data%gru(iGRU)%var(iLookBVAR%routingRunoffFuture)%dat/), start=(/iGRU/),count=(/1,nTimeDelay/)) end do ! iGRU loop diff --git a/build/source/netcdf/read_icond.f90 b/build/source/netcdf/read_icond.f90 index 39267112a..6481856c1 100755 --- a/build/source/netcdf/read_icond.f90 +++ b/build/source/netcdf/read_icond.f90 @@ -21,7 +21,8 @@ module read_icond_module USE nrtype USE netcdf -USE globalData,only:ixHRUfile_min,ixHRUfile_max +USE globalData,only: ixHRUfile_min,ixHRUfile_max +USE globalData,only: nTimeDelay ! number of hours in the time delay histogram implicit none private public::read_icond @@ -169,6 +170,8 @@ subroutine read_icond(iconFile, & ! intent(in): name of USE data_types,only:var_info ! metadata USE get_ixName_module,only:get_varTypeName ! to access type strings for error messages USE updatState_module,only:updateSoil ! update soil states + USE globalData,only:data_step ! timestep of model data / default run + implicit none ! -------------------------------------------------------------------------------------------------------- @@ -188,7 +191,7 @@ subroutine read_icond(iconFile, & ! intent(in): name of integer(i4b) :: fileHRU ! number of HRUs in file integer(i4b) :: fileGRU ! number of GRUs in file integer(i4b) :: iVar, i ! loop indices - integer(i4b),dimension(2) :: ndx ! intermediate array of loop indices + integer(i4b),dimension(1) :: ndx ! intermediate array of loop indices integer(i4b) :: iGRU ! loop index integer(i4b) :: iHRU ! loop index integer(i4b) :: iTDH ! loop index (maybe not needed) @@ -219,9 +222,6 @@ subroutine read_icond(iconFile, & ! intent(in): name of ! Start procedure here err=0; message="read_icond/" - ! set length of time delay histogram - nTDH = 2000 ! TODO figure a way to get this param from allocGlobal into here - ! -------------------------------------------------------------------------------------------------------- ! (1) read the file ! -------------------------------------------------------------------------------------------------------- @@ -240,7 +240,13 @@ subroutine read_icond(iconFile, & ! intent(in): name of ! get dimension of time delay histogram (TDH) in file err = nf90_inq_dimid(ncID,"tdh",dimID); if(err/=nf90_noerr)then; message=trim(message)//'problem finding tdh dimension/'//trim(nf90_strerror(err)); return; end if err = nf90_inquire_dimension(ncID,dimID,len=nTDH); if(err/=nf90_noerr)then; message=trim(message)//'problem reading tdh dimension/'//trim(nf90_strerror(err)); return; end if - ! need check via hardwired value elsewhere in code + + ! check vs hardwired value set in globalData.f90 + if(nTDH /= nTimeDelay)then + write(*,*) 'tdh=',nTDH,' nTimeDelay=',nTimeDelay + message=trim(message)//': state file time delay dimension tdh does not match summa expectation of nTimeDelay set in globalData()' + return + endif ! loop through prognostic variables do iVar = 1,size(prog_meta) @@ -286,7 +292,7 @@ subroutine read_icond(iconFile, & ! intent(in): name of ! get data err = nf90_get_var(ncID,ncVarID,varData); call netcdf_err(err,message) - if(err/=0)then; message=trim(message)//': problem getting the data'; return; endif + if(err/=0)then; message=trim(message)//': problem getting the data for variable '//trim(prog_meta(iVar)%varName); return; endif ! store data in prognostics structure ! loop through GRUs @@ -410,7 +416,7 @@ subroutine read_icond(iconFile, & ! intent(in): name of if(restartFileType/=singleHRU)then ! loop through specific basin variables - ndx = (/iLookBVAR%routingRunoffFuture, iLookBVAR%routingFractionFuture/) ! array of desired variable indices + ndx = (/iLookBVAR%routingRunoffFuture/) ! array of desired variable indices (if more than one eventually) do i = 1,size(ndx) iVar = ndx(i) From c0ff6b5631bc748398b74f587c38b641e637614d Mon Sep 17 00:00:00 2001 From: Andy Wood Date: Wed, 23 Sep 2020 13:33:25 -0600 Subject: [PATCH 04/13] making statefile fix backward compatible --- build/source/driver/summa_restart.f90 | 3 +- build/source/dshare/globalData.f90 | 2 +- build/source/engine/allocspace.f90 | 2 +- build/source/engine/updatState.f90 | 1 - build/source/engine/var_derive.f90 | 4 +- build/source/netcdf/read_icond.f90 | 109 +++++++++++++------------- 6 files changed, 62 insertions(+), 59 deletions(-) diff --git a/build/source/driver/summa_restart.f90 b/build/source/driver/summa_restart.f90 index 5080ec33c..61b80816e 100755 --- a/build/source/driver/summa_restart.f90 +++ b/build/source/driver/summa_restart.f90 @@ -101,13 +101,14 @@ subroutine summa_readRestart(summa1_struc, err, message) nHRU => summa1_struc%nHRU & ! number of global hydrologic response units ) ! assignment to variables in the data structures + ! --------------------------------------------------------------------------------------- ! initialize error control err=0; message='summa_readRestart/' ! identify the start of the writing call date_and_time(values=startRestart) - + ! ***************************************************************************** ! *** read/check initial conditions ! ***************************************************************************** diff --git a/build/source/dshare/globalData.f90 b/build/source/dshare/globalData.f90 index 7ff1a9ce7..efd0c791a 100755 --- a/build/source/dshare/globalData.f90 +++ b/build/source/dshare/globalData.f90 @@ -337,7 +337,7 @@ MODULE globalData integer(i4b),parameter,public :: utcTime=2 ! all times in UTC (timeOffset = longitude/15. hours) integer(i4b),parameter,public :: localTime=3 ! all times local (timeOffset = 0) - ! define fixed dimensions (move these to globalData module?) + ! define fixed dimensions integer(i4b),parameter,public :: nBand=2 ! number of spectral bands integer(i4b),parameter,public :: nTimeDelay=2000 ! number of hours in the time delay histogram (default: ~1 season = 24*365/4) diff --git a/build/source/engine/allocspace.f90 b/build/source/engine/allocspace.f90 index 07a28d643..27e73300a 100755 --- a/build/source/engine/allocspace.f90 +++ b/build/source/engine/allocspace.f90 @@ -550,7 +550,7 @@ subroutine allocateDat_dp(metadata,nSnow,nSoil,nLayers, & ! input ! get the number of variables in the metadata structure nVars = size(metadata) - + ! loop through variables in the data structure do iVar=1,nVars diff --git a/build/source/engine/updatState.f90 b/build/source/engine/updatState.f90 index 700cea109..698c8b1cd 100755 --- a/build/source/engine/updatState.f90 +++ b/build/source/engine/updatState.f90 @@ -70,7 +70,6 @@ subroutine updateSnow(& mLayerVolFracLiq = fLiq*mLayerTheta mLayerVolFracIce = (1._dp - fLiq)*mLayerTheta*(iden_water/iden_ice) !print*, 'mLayerTheta - (mLayerVolFracIce*(iden_ice/iden_water) + mLayerVolFracLiq) = ', mLayerTheta - (mLayerVolFracIce*(iden_ice/iden_water) + mLayerVolFracLiq) - !write(*,'(a,1x,4(f20.10,1x))') 'in updateSnow: fLiq, mLayerTheta, mLayerVolFracIce = ', & ! fLiq, mLayerTheta, mLayerVolFracIce !pause diff --git a/build/source/engine/var_derive.f90 b/build/source/engine/var_derive.f90 index 1b2fbcda2..8227b0407 100755 --- a/build/source/engine/var_derive.f90 +++ b/build/source/engine/var_derive.f90 @@ -418,7 +418,7 @@ subroutine fracFuture(bpar_data,bvar_data,err,message) ! identify number of points in the time-delay runoff variable (should be allocated match nTimeDelay) nTDH = size(runoffFuture) - ! initialize runoffFuture + ! initialize runoffFuture (will be overwritten by initial conditions file values if present) runoffFuture(1:nTDH) = 0._dp ! select option for sub-grid routing @@ -458,7 +458,7 @@ subroutine fracFuture(bpar_data,bvar_data,err,message) if(abs(1._dp - sumFrac) > tolerFrac)then write(*,*) 'fraction of basin runoff histogram being accounted for by time delay vector is ', sumFrac write(*,*) 'this is less than allowed by tolerFrac = ', tolerFrac - message=trim(message)//'not enough bins for the time delay histogram -- fix hard-coded parameter in allocspace.f90' + message=trim(message)//'not enough bins for the time delay histogram -- fix hard-coded parameter in globalData.f90' err=20; return end if ! ensure the fraction sums to one diff --git a/build/source/netcdf/read_icond.f90 b/build/source/netcdf/read_icond.f90 index 6481856c1..443568471 100755 --- a/build/source/netcdf/read_icond.f90 +++ b/build/source/netcdf/read_icond.f90 @@ -194,7 +194,6 @@ subroutine read_icond(iconFile, & ! intent(in): name of integer(i4b),dimension(1) :: ndx ! intermediate array of loop indices integer(i4b) :: iGRU ! loop index integer(i4b) :: iHRU ! loop index - integer(i4b) :: iTDH ! loop index (maybe not needed) integer(i4b) :: dimID ! varible dimension ids integer(i4b) :: ncVarID ! variable ID in netcdf file character(256) :: dimName ! not used except as a placeholder in call to inq_dim function @@ -207,21 +206,20 @@ subroutine read_icond(iconFile, & ! intent(in): name of integer(i4b) :: nSoil, nSnow, nToto ! # layers integer(i4b) :: nTDH ! number of points in time-delay histogram integer(i4b) :: iLayer,jLayer ! layer indices - integer(i4b),parameter :: nBand=2 ! number of spectral bands + integer(i4b),parameter :: nBand=2 ! number of spectral bands integer(i4b) :: nProgVars ! number of prognostic variables written to state file character(len=32),parameter :: scalDimName ='scalarv' ! dimension name for scalar data character(len=32),parameter :: midSoilDimName='midSoil' ! dimension name for soil-only layers character(len=32),parameter :: midTotoDimName='midToto' ! dimension name for layered varaiables character(len=32),parameter :: ifcTotoDimName='ifcToto' ! dimension name for layered varaiables - character(len=32),parameter :: tdhDimName ='tdh' ! dimension name for time-delay basin variables (maybe not needed) + character(len=32),parameter :: tdhDimName ='tdh' ! dimension name for time-delay basin variables - ! -------------------------------------------------------------------------------------------------------- ! Start procedure here err=0; message="read_icond/" - + ! -------------------------------------------------------------------------------------------------------- ! (1) read the file ! -------------------------------------------------------------------------------------------------------- @@ -236,17 +234,6 @@ subroutine read_icond(iconFile, & ! intent(in): name of ! get number of GRUs in file err = nf90_inq_dimid(ncID,"gru",dimID); if(err/=nf90_noerr)then; message=trim(message)//'problem finding gru dimension/'//trim(nf90_strerror(err)); return; end if err = nf90_inquire_dimension(ncID,dimID,len=fileGRU); if(err/=nf90_noerr)then; message=trim(message)//'problem reading gru dimension/'//trim(nf90_strerror(err)); return; end if - - ! get dimension of time delay histogram (TDH) in file - err = nf90_inq_dimid(ncID,"tdh",dimID); if(err/=nf90_noerr)then; message=trim(message)//'problem finding tdh dimension/'//trim(nf90_strerror(err)); return; end if - err = nf90_inquire_dimension(ncID,dimID,len=nTDH); if(err/=nf90_noerr)then; message=trim(message)//'problem reading tdh dimension/'//trim(nf90_strerror(err)); return; end if - - ! check vs hardwired value set in globalData.f90 - if(nTDH /= nTimeDelay)then - write(*,*) 'tdh=',nTDH,' nTimeDelay=',nTimeDelay - message=trim(message)//': state file time delay dimension tdh does not match summa expectation of nTimeDelay set in globalData()' - return - endif ! loop through prognostic variables do iVar = 1,size(prog_meta) @@ -409,57 +396,73 @@ subroutine read_icond(iconFile, & ! intent(in): name of end do ! looping through GRUs ! -------------------------------------------------------------------------------------------------------- - ! (2) now get the two basin runoff variables + ! (2) now get the basin variable(s) ! -------------------------------------------------------------------------------------------------------- ! get the index in the file: single HRU if(restartFileType/=singleHRU)then + + ! get dimension of time delay histogram (TDH) from initial conditions file + err = nf90_inq_dimid(ncID,"tdh",dimID); + if(err/=nf90_noerr)then + write(*,*) 'WARNING: routingRunoffFuture is not in the initial conditions file ... using zeros' ! previously created in var_derive.f90 + err=nf90_noerr ! reset this err + + else + err = nf90_inquire_dimension(ncID,dimID,len=nTDH); + if(err/=nf90_noerr)then; message=trim(message)//'problem reading tdh dimension from initial condition file/'//trim(nf90_strerror(err)); return; end if + + ! check vs hardwired value set in globalData.f90 + if(nTDH /= nTimeDelay)then + write(*,*) 'tdh=',nTDH,' nTimeDelay=',nTimeDelay + message=trim(message)//': state file time delay dimension tdh does not match summa expectation of nTimeDelay set in globalData()' + return + endif - ! loop through specific basin variables - ndx = (/iLookBVAR%routingRunoffFuture/) ! array of desired variable indices (if more than one eventually) - do i = 1,size(ndx) - iVar = ndx(i) + ! loop through specific basin variables (currently 1 but loop provided to enable inclusion of others) + ndx = (/iLookBVAR%routingRunoffFuture/) ! array of desired variable indices + do i = 1,size(ndx) + iVar = ndx(i) - ! get tdh dimension Id in file (should be 'tdh') - err = nf90_inq_dimid(ncID,trim(tdhDimName), dimID); - if(err/=0)then; message=trim(message)//': problem with dimension ids for tdh vars'; return; endif + ! get tdh dimension Id in file (should be 'tdh') + err = nf90_inq_dimid(ncID,trim(tdhDimName), dimID); + if(err/=0)then; message=trim(message)//': problem with dimension ids for tdh vars'; return; endif - ! get the tdh dimension length (dimName and dimLen are outputs of this call) - err = nf90_inquire_dimension(ncID,dimID,dimName,dimLen); call netcdf_err(err,message) - if(err/=0)then; message=trim(message)//': problem getting the dimension length for tdh vars'; return; endif + ! get the tdh dimension length (dimName and dimLen are outputs of this call) + err = nf90_inquire_dimension(ncID,dimID,dimName,dimLen); call netcdf_err(err,message) + if(err/=0)then; message=trim(message)//': problem getting the dimension length for tdh vars'; return; endif - ! get tdh-based variable id - err = nf90_inq_varid(ncID,trim(bvar_meta(iVar)%varName),ncVarID); call netcdf_err(err,message) - if(err/=0)then; message=trim(message)//': problem with getting basin variable id, var='//trim(bvar_meta(iVar)%varName); return; endif + ! get tdh-based variable id + err = nf90_inq_varid(ncID,trim(bvar_meta(iVar)%varName),ncVarID); call netcdf_err(err,message) + if(err/=0)then; message=trim(message)//': problem with getting basin variable id, var='//trim(bvar_meta(iVar)%varName); return; endif - ! initialize the tdh variable data - allocate(varData(fileGRU,dimLen),stat=err) - if(err/=0)then; print*, 'err= ',err; message=trim(message)//'problem allocating GRU variable data'; return; endif - - ! get data - err = nf90_get_var(ncID,ncVarID,varData); call netcdf_err(err,message) - if(err/=0)then; message=trim(message)//': problem getting the data'; return; endif - - ! store data in basin var (bvar) structure - do iGRU = 1,nGRU - - ! put the data into data structures - bvarData%gru(iGRU)%var(iVar)%dat(1:nTDH) = varData(iGRU,1:nTDH) - ! check whether the first values is set to nf90_fill_double - if(any(abs(bvarData%gru(iGRU)%var(iVar)%dat(1:nTDH) - nf90_fill_double) < epsilon(varData)))then; err=20; endif - if(err==20)then; message=trim(message)//"data set to the fill value (name='"//trim(bvar_meta(iVar)%varName)//"')"; return; endif + ! initialize the tdh variable data + allocate(varData(fileGRU,dimLen),stat=err) + if(err/=0)then; print*, 'err= ',err; message=trim(message)//'problem allocating GRU variable data'; return; endif + + ! get data + err = nf90_get_var(ncID,ncVarID,varData); call netcdf_err(err,message) + if(err/=0)then; message=trim(message)//': problem getting the data'; return; endif + + ! store data in basin var (bvar) structure + do iGRU = 1,nGRU + + ! put the data into data structures + bvarData%gru(iGRU)%var(iVar)%dat(1:nTDH) = varData(iGRU,1:nTDH) + ! check whether the first values is set to nf90_fill_double + if(any(abs(bvarData%gru(iGRU)%var(iVar)%dat(1:nTDH) - nf90_fill_double) < epsilon(varData)))then; err=20; endif + if(err==20)then; message=trim(message)//"data set to the fill value (name='"//trim(bvar_meta(iVar)%varName)//"')"; return; endif - end do ! end iGRU loop + end do ! end iGRU loop - ! deallocate temporary data array for next variable - deallocate(varData, stat=err) - if(err/=0)then; message=trim(message)//'problem deallocating GRU variable data'; return; endif + ! deallocate temporary data array for next variable + deallocate(varData, stat=err) + if(err/=0)then; message=trim(message)//'problem deallocating GRU variable data'; return; endif - end do ! end looping through basin variables - + end do ! end looping through basin variables + endif ! end if case for tdh variables being in init. cond. file endif ! end if case for not being a singleHRU run - end subroutine read_icond end module read_icond_module From c5b5240d14f07de90f92bc83e7bca635ae201fee Mon Sep 17 00:00:00 2001 From: Andy Wood Date: Wed, 23 Sep 2020 16:42:01 -0600 Subject: [PATCH 05/13] further cleanup on statefile fix --- build/source/driver/summa_writeOutput.f90 | 1 - build/source/netcdf/read_icond.f90 | 4 +--- 2 files changed, 1 insertion(+), 4 deletions(-) mode change 100755 => 100644 build/source/netcdf/read_icond.f90 diff --git a/build/source/driver/summa_writeOutput.f90 b/build/source/driver/summa_writeOutput.f90 index 121fde7bc..32ef3d466 100755 --- a/build/source/driver/summa_writeOutput.f90 +++ b/build/source/driver/summa_writeOutput.f90 @@ -287,7 +287,6 @@ subroutine summa_writeOutputFiles(modelTimeStep, summa1_struc, err, message) restartFile=trim(STATE_PATH)//trim(OUTPUT_PREFIX)//'_restart_'//trim(timeString)//trim(output_fileSuffix)//'.nc' endif - !call writeRestart(restartFile,nGRU,nHRU,prog_meta,progStruct,maxLayers,maxSnowLayers,indx_meta,indxStruct,err,cmessage) call writeRestart(restartFile,nGRU,nHRU,prog_meta,progStruct,bvar_meta,bvarStruct,maxLayers,maxSnowLayers,indx_meta,indxStruct,err,cmessage) if(err/=0)then; message=trim(message)//trim(cmessage); return; endif end if diff --git a/build/source/netcdf/read_icond.f90 b/build/source/netcdf/read_icond.f90 old mode 100755 new mode 100644 index 443568471..29b4859b0 --- a/build/source/netcdf/read_icond.f90 +++ b/build/source/netcdf/read_icond.f90 @@ -170,7 +170,6 @@ subroutine read_icond(iconFile, & ! intent(in): name of USE data_types,only:var_info ! metadata USE get_ixName_module,only:get_varTypeName ! to access type strings for error messages USE updatState_module,only:updateSoil ! update soil states - USE globalData,only:data_step ! timestep of model data / default run implicit none @@ -391,8 +390,7 @@ subroutine read_icond(iconFile, & ! intent(in): name of if (err/=0) then; message=trim(message)//trim(cmessage); return; end if end do ! looping through soil layers - - end do ! looping throuygh HRUs + end do ! looping through HRUs end do ! looping through GRUs ! -------------------------------------------------------------------------------------------------------- From 462509b13fbfb2b8e1e04e17cef40e62d1a38847 Mon Sep 17 00:00:00 2001 From: Andy Wood Date: Wed, 23 Sep 2020 17:47:06 -0600 Subject: [PATCH 06/13] correcting check for gru dim in state file if tdh vars not present --- build/source/netcdf/read_icond.f90 | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/build/source/netcdf/read_icond.f90 b/build/source/netcdf/read_icond.f90 index 29b4859b0..ec84a9ffb 100644 --- a/build/source/netcdf/read_icond.f90 +++ b/build/source/netcdf/read_icond.f90 @@ -229,10 +229,6 @@ subroutine read_icond(iconFile, & ! intent(in): name of ! get number of HRUs in file err = nf90_inq_dimid(ncID,"hru",dimID); if(err/=nf90_noerr)then; message=trim(message)//'problem finding hru dimension/'//trim(nf90_strerror(err)); return; end if err = nf90_inquire_dimension(ncID,dimID,len=fileHRU); if(err/=nf90_noerr)then; message=trim(message)//'problem reading hru dimension/'//trim(nf90_strerror(err)); return; end if - - ! get number of GRUs in file - err = nf90_inq_dimid(ncID,"gru",dimID); if(err/=nf90_noerr)then; message=trim(message)//'problem finding gru dimension/'//trim(nf90_strerror(err)); return; end if - err = nf90_inquire_dimension(ncID,dimID,len=fileGRU); if(err/=nf90_noerr)then; message=trim(message)//'problem reading gru dimension/'//trim(nf90_strerror(err)); return; end if ! loop through prognostic variables do iVar = 1,size(prog_meta) @@ -405,11 +401,16 @@ subroutine read_icond(iconFile, & ! intent(in): name of if(err/=nf90_noerr)then write(*,*) 'WARNING: routingRunoffFuture is not in the initial conditions file ... using zeros' ! previously created in var_derive.f90 err=nf90_noerr ! reset this err - + else + ! the state file *does* have the basin variable(s), so process them err = nf90_inquire_dimension(ncID,dimID,len=nTDH); if(err/=nf90_noerr)then; message=trim(message)//'problem reading tdh dimension from initial condition file/'//trim(nf90_strerror(err)); return; end if + ! get number of GRUs in file + err = nf90_inq_dimid(ncID,"gru",dimID); if(err/=nf90_noerr)then; message=trim(message)//'problem finding gru dimension/'//trim(nf90_strerror(err)); return; end if + err = nf90_inquire_dimension(ncID,dimID,len=fileGRU); if(err/=nf90_noerr)then; message=trim(message)//'problem reading gru dimension/'//trim(nf90_strerror(err)); return; end if + ! check vs hardwired value set in globalData.f90 if(nTDH /= nTimeDelay)then write(*,*) 'tdh=',nTDH,' nTimeDelay=',nTimeDelay From fb2382092507c06a818e96d5f0240cdfa56922e3 Mon Sep 17 00:00:00 2001 From: arbennett Date: Thu, 24 Sep 2020 11:55:40 -0700 Subject: [PATCH 07/13] Add potential fix? --- build/source/netcdf/read_icond.f90 | 36 +++++++++++++++++------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/build/source/netcdf/read_icond.f90 b/build/source/netcdf/read_icond.f90 index ec84a9ffb..fdf491c9c 100644 --- a/build/source/netcdf/read_icond.f90 +++ b/build/source/netcdf/read_icond.f90 @@ -42,6 +42,7 @@ subroutine read_icond_nlayers(iconFile,nGRU,indx_meta,err,message) USE nrtype USE var_lookup,only:iLookIndex ! variable lookup structure USE globalData,only:gru_struc ! gru-hru mapping structures + USE globalData,only:startGRU ! index of first gru for parallel runs USE netcdf_util_module,only:nc_file_close ! close netcdf file USE netcdf_util_module,only:nc_file_open ! close netcdf file USE netcdf_util_module,only:netcdf_err ! netcdf error handling @@ -159,6 +160,7 @@ subroutine read_icond(iconFile, & ! intent(in): name of USE globalData,only:prog_meta ! metadata for prognostic variables USE globalData,only:bvar_meta ! metadata for basin (GRU) variables USE globalData,only:gru_struc ! gru-hru mapping structures + USE globalData,only:startGRU ! index of first gru for parallel runs USE globaldata,only:iname_soil,iname_snow ! named variables to describe the type of layer USE netcdf_util_module,only:nc_file_open ! open netcdf file USE netcdf_util_module,only:nc_file_close ! close netcdf file @@ -213,7 +215,7 @@ subroutine read_icond(iconFile, & ! intent(in): name of character(len=32),parameter :: midTotoDimName='midToto' ! dimension name for layered varaiables character(len=32),parameter :: ifcTotoDimName='ifcToto' ! dimension name for layered varaiables character(len=32),parameter :: tdhDimName ='tdh' ! dimension name for time-delay basin variables - + ! -------------------------------------------------------------------------------------------------------- ! Start procedure here @@ -229,7 +231,7 @@ subroutine read_icond(iconFile, & ! intent(in): name of ! get number of HRUs in file err = nf90_inq_dimid(ncID,"hru",dimID); if(err/=nf90_noerr)then; message=trim(message)//'problem finding hru dimension/'//trim(nf90_strerror(err)); return; end if err = nf90_inquire_dimension(ncID,dimID,len=fileHRU); if(err/=nf90_noerr)then; message=trim(message)//'problem reading hru dimension/'//trim(nf90_strerror(err)); return; end if - + ! loop through prognostic variables do iVar = 1,size(prog_meta) @@ -274,6 +276,7 @@ subroutine read_icond(iconFile, & ! intent(in): name of ! get data err = nf90_get_var(ncID,ncVarID,varData); call netcdf_err(err,message) + write(*,*) shape(varData), shape(gru_struc), ixFile if(err/=0)then; message=trim(message)//': problem getting the data for variable '//trim(prog_meta(iVar)%varName); return; endif ! store data in prognostics structure @@ -283,7 +286,7 @@ subroutine read_icond(iconFile, & ! intent(in): name of iHRU_global = gru_struc(iGRU)%hruInfo(iHRU)%hru_nc iHRU_local = (iHRU_global - ixHRUfile_min) + 1 - + ! get the number of layers nSnow = gru_struc(iGRU)%hruInfo(iHRU)%nSnow nSoil = gru_struc(iGRU)%hruInfo(iHRU)%nSoil @@ -291,11 +294,12 @@ subroutine read_icond(iconFile, & ! intent(in): name of ! get the index in the file: single HRU if(restartFileType==singleHRU)then - ixFile = 1 ! use for single HRU restart file + ! ixFile = 1 ! use for single HRU restart file + ixFile = startGRU ! get the index in the file: multi HRU else - ixFile = iHRU_local + ixFile = startGRU + iHRU_local - 1 endif ! put the data into data structures and check that none of the values are set to nf90_fill_double @@ -335,7 +339,7 @@ subroutine read_icond(iconFile, & ! intent(in): name of if(err/=0)then; message=trim(message)//'problem deallocating HRU variable data'; return; endif end do ! end looping through prognostic variables (iVar) - + ! -------------------------------------------------------------------------------------------------------- ! (2) set number of layers ! -------------------------------------------------------------------------------------------------------- @@ -388,20 +392,20 @@ subroutine read_icond(iconFile, & ! intent(in): name of end do ! looping through soil layers end do ! looping through HRUs end do ! looping through GRUs - + ! -------------------------------------------------------------------------------------------------------- ! (2) now get the basin variable(s) ! -------------------------------------------------------------------------------------------------------- ! get the index in the file: single HRU if(restartFileType/=singleHRU)then - + ! get dimension of time delay histogram (TDH) from initial conditions file err = nf90_inq_dimid(ncID,"tdh",dimID); if(err/=nf90_noerr)then write(*,*) 'WARNING: routingRunoffFuture is not in the initial conditions file ... using zeros' ! previously created in var_derive.f90 err=nf90_noerr ! reset this err - + else ! the state file *does* have the basin variable(s), so process them err = nf90_inquire_dimension(ncID,dimID,len=nTDH); @@ -422,19 +426,19 @@ subroutine read_icond(iconFile, & ! intent(in): name of ndx = (/iLookBVAR%routingRunoffFuture/) ! array of desired variable indices do i = 1,size(ndx) iVar = ndx(i) - + ! get tdh dimension Id in file (should be 'tdh') - err = nf90_inq_dimid(ncID,trim(tdhDimName), dimID); + err = nf90_inq_dimid(ncID,trim(tdhDimName), dimID); if(err/=0)then; message=trim(message)//': problem with dimension ids for tdh vars'; return; endif ! get the tdh dimension length (dimName and dimLen are outputs of this call) err = nf90_inquire_dimension(ncID,dimID,dimName,dimLen); call netcdf_err(err,message) if(err/=0)then; message=trim(message)//': problem getting the dimension length for tdh vars'; return; endif - + ! get tdh-based variable id err = nf90_inq_varid(ncID,trim(bvar_meta(iVar)%varName),ncVarID); call netcdf_err(err,message) if(err/=0)then; message=trim(message)//': problem with getting basin variable id, var='//trim(bvar_meta(iVar)%varName); return; endif - + ! initialize the tdh variable data allocate(varData(fileGRU,dimLen),stat=err) if(err/=0)then; print*, 'err= ',err; message=trim(message)//'problem allocating GRU variable data'; return; endif @@ -451,13 +455,13 @@ subroutine read_icond(iconFile, & ! intent(in): name of ! check whether the first values is set to nf90_fill_double if(any(abs(bvarData%gru(iGRU)%var(iVar)%dat(1:nTDH) - nf90_fill_double) < epsilon(varData)))then; err=20; endif if(err==20)then; message=trim(message)//"data set to the fill value (name='"//trim(bvar_meta(iVar)%varName)//"')"; return; endif - + end do ! end iGRU loop - + ! deallocate temporary data array for next variable deallocate(varData, stat=err) if(err/=0)then; message=trim(message)//'problem deallocating GRU variable data'; return; endif - + end do ! end looping through basin variables endif ! end if case for tdh variables being in init. cond. file endif ! end if case for not being a singleHRU run From ee5eff21c73799deb4b0f327b8171a91950186c4 Mon Sep 17 00:00:00 2001 From: Andy Wood Date: Thu, 24 Sep 2020 13:12:40 -0600 Subject: [PATCH 08/13] misc state file fix work in progress --- build/source/dshare/globalData.f90 | 2 +- build/source/netcdf/read_icond.f90 | 23 +++++++++++++---------- 2 files changed, 14 insertions(+), 11 deletions(-) diff --git a/build/source/dshare/globalData.f90 b/build/source/dshare/globalData.f90 index efd0c791a..68300c427 100755 --- a/build/source/dshare/globalData.f90 +++ b/build/source/dshare/globalData.f90 @@ -250,7 +250,7 @@ MODULE globalData ! define metadata for model forcing datafile type(file_info),save,public,allocatable :: forcFileInfo(:) ! file info for model forcing data - ! define indices describing the indices of the first and last HRUs in the forcing file + ! define index variables describing the indices of the first and last HRUs in the forcing file integer(i4b),save,public :: ixHRUfile_min ! minimum index integer(i4b),save,public :: ixHRUfile_max ! maximum index diff --git a/build/source/netcdf/read_icond.f90 b/build/source/netcdf/read_icond.f90 index ec84a9ffb..aa1fb662c 100644 --- a/build/source/netcdf/read_icond.f90 +++ b/build/source/netcdf/read_icond.f90 @@ -79,29 +79,27 @@ subroutine read_icond_nlayers(iconFile,nGRU,indx_meta,err,message) call nc_file_open(iconFile,nf90_nowrite,ncid,err,cmessage); if (err/=0) then; message=trim(message)//trim(cmessage); return; end if - ! get number of HRUs in file + ! get number of HRUs in file (the GRU variable(s), if present, are processed at the end) err = nf90_inq_dimid(ncID,"hru",dimId); if(err/=nf90_noerr)then; message=trim(message)//'problem finding hru dimension/'//trim(nf90_strerror(err)); return; end if err = nf90_inquire_dimension(ncID,dimId,len=fileHRU); if(err/=nf90_noerr)then; message=trim(message)//'problem reading hru dimension/'//trim(nf90_strerror(err)); return; end if - ! allocate storage for reading from file + ! allocate storage for reading from file (allocate entire file size, even when doing subdomain run) allocate(snowData(fileHRU)) allocate(soilData(fileHRU)) snowData = 0 soilData = 0 - ! get variable ids + ! get netcdf ids for the variables holding number of snow and soil layers in each hru err = nf90_inq_varid(ncid,trim(indx_meta(iLookIndex%nSnow)%varName),snowid); call netcdf_err(err,message) err = nf90_inq_varid(ncid,trim(indx_meta(iLookIndex%nSoil)%varName),soilid); call netcdf_err(err,message) - ! get data + ! get nSnow and nSoil data (reads entire state file) err = nf90_get_var(ncid,snowid,snowData); call netcdf_err(err,message) err = nf90_get_var(ncid,soilid,soilData); call netcdf_err(err,message) - !print*, 'snowData = ', snowData - !print*, 'soilData = ', soilData ixHRUfile_min=huge(1) ixHRUfile_max=0 - ! assign to index structure - gru by hru + ! find the min and max hru indices in the state file (why? this should be known from the ) do iGRU = 1,nGRU do iHRU = 1,gru_struc(iGRU)%hruCount if(gru_struc(iGRU)%hruInfo(iHRU)%hru_nc < ixHRUfile_min) ixHRUfile_min = gru_struc(iGRU)%hruInfo(iHRU)%hru_nc @@ -109,20 +107,25 @@ subroutine read_icond_nlayers(iconFile,nGRU,indx_meta,err,message) end do end do + ! tmp-note: nGRU = size(gru_struc) ... has number of grus in requested gru subset run + ! iHRU_global contains the correct subset index do iGRU = 1,nGRU do iHRU = 1,gru_struc(iGRU)%hruCount iHRU_global = gru_struc(iGRU)%hruInfo(iHRU)%hru_nc iHRU_local = (iHRU_global - ixHRUfile_min) + 1 + print*, 'iGRU, HRU, iHRU_global, iHRU_local = ', iGRU, iHRU, iHRU_global, iHRU_local ! single HRU - if(restartFileType==singleHRU)then + if(restartFileType==singleHRU)then ! restartFileType hardwired above (to multiHRU) gru_struc(iGRU)%hruInfo(iHRU)%nSnow = snowData(1) gru_struc(iGRU)%hruInfo(iHRU)%nSoil = soilData(1) ! multi HRU else - gru_struc(iGRU)%hruInfo(iHRU)%nSnow = snowData(iHRU_local) - gru_struc(iGRU)%hruInfo(iHRU)%nSoil = soilData(iHRU_local) + !gru_struc(iGRU)%hruInfo(iHRU)%nSnow = snowData(iHRU_local) + !gru_struc(iGRU)%hruInfo(iHRU)%nSoil = soilData(iHRU_local) + gru_struc(iGRU)%hruInfo(iHRU)%nSnow = snowData(iHRU_global) + gru_struc(iGRU)%hruInfo(iHRU)%nSoil = soilData(iHRU_global) endif end do From c7a19af2ee7f6b5644aa13d13ea9573ade1cd2b3 Mon Sep 17 00:00:00 2001 From: arbennett Date: Thu, 24 Sep 2020 13:01:33 -0700 Subject: [PATCH 09/13] Fix restarts on basin vars when using -g --- build/source/netcdf/read_icond.f90 | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/build/source/netcdf/read_icond.f90 b/build/source/netcdf/read_icond.f90 index fdf491c9c..753c538c9 100644 --- a/build/source/netcdf/read_icond.f90 +++ b/build/source/netcdf/read_icond.f90 @@ -276,7 +276,6 @@ subroutine read_icond(iconFile, & ! intent(in): name of ! get data err = nf90_get_var(ncID,ncVarID,varData); call netcdf_err(err,message) - write(*,*) shape(varData), shape(gru_struc), ixFile if(err/=0)then; message=trim(message)//': problem getting the data for variable '//trim(prog_meta(iVar)%varName); return; endif ! store data in prognostics structure @@ -294,9 +293,7 @@ subroutine read_icond(iconFile, & ! intent(in): name of ! get the index in the file: single HRU if(restartFileType==singleHRU)then - ! ixFile = 1 ! use for single HRU restart file - ixFile = startGRU - + ixFile = 1 ! use for single HRU restart file ! get the index in the file: multi HRU else ixFile = startGRU + iHRU_local - 1 @@ -451,7 +448,7 @@ subroutine read_icond(iconFile, & ! intent(in): name of do iGRU = 1,nGRU ! put the data into data structures - bvarData%gru(iGRU)%var(iVar)%dat(1:nTDH) = varData(iGRU,1:nTDH) + bvarData%gru(iGRU)%var(iVar)%dat(1:nTDH) = varData((iGRU+startGRU-1),1:nTDH) ! check whether the first values is set to nf90_fill_double if(any(abs(bvarData%gru(iGRU)%var(iVar)%dat(1:nTDH) - nf90_fill_double) < epsilon(varData)))then; err=20; endif if(err==20)then; message=trim(message)//"data set to the fill value (name='"//trim(bvar_meta(iVar)%varName)//"')"; return; endif From 72a1019193728dc2754473b75a3b2c6a77fefd96 Mon Sep 17 00:00:00 2001 From: Andy Wood Date: Thu, 24 Sep 2020 14:35:46 -0600 Subject: [PATCH 10/13] remove debug print statement in read_icond.f90 --- build/source/netcdf/read_icond.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/build/source/netcdf/read_icond.f90 b/build/source/netcdf/read_icond.f90 index 4d21349b0..b5a30632a 100644 --- a/build/source/netcdf/read_icond.f90 +++ b/build/source/netcdf/read_icond.f90 @@ -114,7 +114,6 @@ subroutine read_icond_nlayers(iconFile,nGRU,indx_meta,err,message) do iHRU = 1,gru_struc(iGRU)%hruCount iHRU_global = gru_struc(iGRU)%hruInfo(iHRU)%hru_nc iHRU_local = (iHRU_global - ixHRUfile_min) + 1 - print*, 'iGRU, HRU, iHRU_global, iHRU_local = ', iGRU, iHRU, iHRU_global, iHRU_local ! single HRU if(restartFileType==singleHRU)then ! restartFileType hardwired above (to multiHRU) From fc7f5f3f6c0618a12d594bcf08acc8b416b51fac Mon Sep 17 00:00:00 2001 From: arbennett Date: Thu, 24 Sep 2020 15:05:46 -0700 Subject: [PATCH 11/13] Output gruid and hruid as coordinates of gru and hru --- build/source/netcdf/def_output.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/build/source/netcdf/def_output.f90 b/build/source/netcdf/def_output.f90 index 792196735..d376012c9 100755 --- a/build/source/netcdf/def_output.f90 +++ b/build/source/netcdf/def_output.f90 @@ -469,14 +469,14 @@ subroutine write_hru_info(ncid, err, message) do iGRU = 1, size(gru_struc) ! GRU info - err = nf90_put_var(ncid, gruVarID, gru_struc(iGRU)%gru_nc, start=(/iGRU/)) + err = nf90_put_var(ncid, gruVarID, gru_struc(iGRU)%gruId, start=(/iGRU/)) if (err/=nf90_NoErr) then; message=trim(message)//'nf90_write_gruVar'; call netcdf_err(err,message); return; end if err = nf90_put_var(ncid, gruIdVarID, gru_struc(iGRU)%gruId, start=(/iGRU/)) if (err/=nf90_NoErr) then; message=trim(message)//'nf90_write_gruIdVar'; call netcdf_err(err,message); return; end if ! HRU info do iHRU = 1, gru_struc(iGRU)%hruCount - err = nf90_put_var(ncid, hruVarID, gru_struc(iGRU)%hruInfo(iHRU)%hru_nc, start=(/gru_struc(iGRU)%hruInfo(iHRU)%hru_ix/)) + err = nf90_put_var(ncid, hruVarID, gru_struc(iGRU)%hruInfo(iHRU)%hru_id, start=(/gru_struc(iGRU)%hruInfo(iHRU)%hru_ix/)) if (err/=nf90_NoErr) then; message=trim(message)//'nf90_write_hruVar'; call netcdf_err(err,message); return; end if err = nf90_put_var(ncid, hruIdVarID, gru_struc(iGRU)%hruInfo(iHRU)%hru_id, start=(/gru_struc(iGRU)%hruInfo(iHRU)%hru_ix/)) if (err/=nf90_NoErr) then; message=trim(message)//'nf90_write_hruIdVar'; call netcdf_err(err,message); return; end if From 5039e9a87f07ad9a991ed886341c9f9bb7707003 Mon Sep 17 00:00:00 2001 From: Andy Wood Date: Wed, 7 Oct 2020 11:48:14 -0600 Subject: [PATCH 12/13] Fixes small bug from PR #433 -- int64 ids were being written to int-defined variables while writing coordinates; Makes variable naming consistent in gru/hru structures (gru_id ~ hru_id); expands gitignore to include pynb tracking files --- build/source/dshare/data_types.f90 | 2 +- build/source/engine/read_attrb.f90 | 6 +++--- build/source/netcdf/def_output.f90 | 12 ++++++------ 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/build/source/dshare/data_types.f90 b/build/source/dshare/data_types.f90 index 086b749ab..cf20b1e89 100755 --- a/build/source/dshare/data_types.f90 +++ b/build/source/dshare/data_types.f90 @@ -112,7 +112,7 @@ MODULE data_types ! define mapping from GRUs to the HRUs type, public :: gru2hru_map - integer(8) :: gruId ! id of the gru + integer(8) :: gru_id ! id of the gru integer(i4b) :: hruCount ! total number of hrus in the gru type(hru_info), allocatable :: hruInfo(:) ! basic information of HRUs within the gru integer(i4b) :: gru_nc ! index of gru in the netcdf file diff --git a/build/source/engine/read_attrb.f90 b/build/source/engine/read_attrb.f90 index 10124a414..e07fff1a1 100755 --- a/build/source/engine/read_attrb.f90 +++ b/build/source/engine/read_attrb.f90 @@ -143,7 +143,7 @@ subroutine read_dimension(attrFile,fileGRU,fileHRU,nGRU,nHRU,err,message,startGR ! gru to hru mapping iGRU = 1 gru_struc(iGRU)%hruCount = 1 ! number of HRUs in each GRU - gru_struc(iGRU)%gruId = hru2gru_id(checkHRU) ! set gru id + gru_struc(iGRU)%gru_id = hru2gru_id(checkHRU) ! set gru id gru_struc(iGRU)%gru_nc = sGRU ! set gru index within the netcdf file allocate(gru_struc(iGRU)%hruInfo(gru_struc(iGRU)%hruCount)) ! allocate second level of gru to hru map gru_struc(iGRU)%hruInfo(iGRU)%hru_nc = checkHRU ! set hru id in attributes netcdf file @@ -157,11 +157,11 @@ subroutine read_dimension(attrFile,fileGRU,fileHRU,nGRU,nHRU,err,message,startGR if (count(hru2gru_Id == gru_id(iGRU+sGRU-1)) < 1) then; err=20; message=trim(message)//'problem finding HRUs belonging to GRU'; return; end if gru_struc(iGRU)%hruCount = count(hru2gru_Id == gru_id(iGRU+sGRU-1)) ! number of HRUs in each GRU - gru_struc(iGRU)%gruId = gru_id(iGRU+sGRU-1) ! set gru id + gru_struc(iGRU)%gru_id = gru_id(iGRU+sGRU-1) ! set gru id gru_struc(iGRU)%gru_nc = iGRU+sGRU-1 ! set gru index in the netcdf file allocate(gru_struc(iGRU)%hruInfo(gru_struc(iGRU)%hruCount)) ! allocate second level of gru to hru map - gru_struc(iGRU)%hruInfo(:)%hru_nc = pack(hru_ix,hru2gru_id == gru_struc(iGRU)%gruId) ! set hru id in attributes netcdf file + gru_struc(iGRU)%hruInfo(:)%hru_nc = pack(hru_ix,hru2gru_id == gru_struc(iGRU)%gru_id) ! set hru id in attributes netcdf file gru_struc(iGRU)%hruInfo(:)%hru_ix = arth(iHRU,1,gru_struc(iGRU)%hruCount) ! set index of hru in run domain gru_struc(iGRU)%hruInfo(:)%hru_id = hru_id(gru_struc(iGRU)%hruInfo(:)%hru_nc) ! set id of hru iHRU = iHRU + gru_struc(iGRU)%hruCount diff --git a/build/source/netcdf/def_output.f90 b/build/source/netcdf/def_output.f90 index d376012c9..996eef1c1 100755 --- a/build/source/netcdf/def_output.f90 +++ b/build/source/netcdf/def_output.f90 @@ -443,13 +443,13 @@ subroutine write_hru_info(ncid, err, message) err = nf90_redef(ncid); call netcdf_err(err, message); if (err/=nf90_NoErr) return ! define HRU var - err = nf90_def_var(ncid, trim(hru_DimName), nf90_int, hru_DimID, hruVarID); if (err/=nf90_NoErr) then; message=trim(message)//'nf90_define_hruVar' ; call netcdf_err(err,message); return; end if - err = nf90_put_att(ncid, hruVarID, 'long_name', 'hru index in the input file'); if (err/=nf90_NoErr) then; message=trim(message)//'write_hruVar_longname'; call netcdf_err(err,message); return; end if + err = nf90_def_var(ncid, trim(hru_DimName), nf90_int64, hru_DimID, hruVarID); if (err/=nf90_NoErr) then; message=trim(message)//'nf90_define_hruVar' ; call netcdf_err(err,message); return; end if + err = nf90_put_att(ncid, hruVarID, 'long_name', 'hruId in the input file'); if (err/=nf90_NoErr) then; message=trim(message)//'write_hruVar_longname'; call netcdf_err(err,message); return; end if err = nf90_put_att(ncid, hruVarID, 'units', '-' ); if (err/=nf90_NoErr) then; message=trim(message)//'write_hruVar_unit'; call netcdf_err(err,message); return; end if ! define GRU var - err = nf90_def_var(ncid, trim(gru_DimName), nf90_int, gru_DimID, gruVarID); if (err/=nf90_NoErr) then; message=trim(message)//'nf90_define_gruVar' ; call netcdf_err(err,message); return; end if - err = nf90_put_att(ncid, gruVarID, 'long_name', 'gru index in the input file'); if (err/=nf90_NoErr) then; message=trim(message)//'write_gruVar_longname'; call netcdf_err(err,message); return; end if + err = nf90_def_var(ncid, trim(gru_DimName), nf90_int64, gru_DimID, gruVarID); if (err/=nf90_NoErr) then; message=trim(message)//'nf90_define_gruVar' ; call netcdf_err(err,message); return; end if + err = nf90_put_att(ncid, gruVarID, 'long_name', 'gruId in the input file'); if (err/=nf90_NoErr) then; message=trim(message)//'write_gruVar_longname'; call netcdf_err(err,message); return; end if err = nf90_put_att(ncid, gruVarID, 'units', '-' ); if (err/=nf90_NoErr) then; message=trim(message)//'write_gruVar_unit'; call netcdf_err(err,message); return; end if ! define hruId var @@ -469,9 +469,9 @@ subroutine write_hru_info(ncid, err, message) do iGRU = 1, size(gru_struc) ! GRU info - err = nf90_put_var(ncid, gruVarID, gru_struc(iGRU)%gruId, start=(/iGRU/)) + err = nf90_put_var(ncid, gruVarID, gru_struc(iGRU)%gru_id, start=(/iGRU/)) if (err/=nf90_NoErr) then; message=trim(message)//'nf90_write_gruVar'; call netcdf_err(err,message); return; end if - err = nf90_put_var(ncid, gruIdVarID, gru_struc(iGRU)%gruId, start=(/iGRU/)) + err = nf90_put_var(ncid, gruIdVarID, gru_struc(iGRU)%gru_id, start=(/iGRU/)) if (err/=nf90_NoErr) then; message=trim(message)//'nf90_write_gruIdVar'; call netcdf_err(err,message); return; end if ! HRU info From e66105f459d69135457211ce03801c5d5f554dfd Mon Sep 17 00:00:00 2001 From: Andy Wood Date: Wed, 7 Oct 2020 12:25:37 -0600 Subject: [PATCH 13/13] cleaning up comments in read_icond.f90 --- build/source/netcdf/read_icond.f90 | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/build/source/netcdf/read_icond.f90 b/build/source/netcdf/read_icond.f90 index b5a30632a..c5bdd929e 100644 --- a/build/source/netcdf/read_icond.f90 +++ b/build/source/netcdf/read_icond.f90 @@ -100,7 +100,7 @@ subroutine read_icond_nlayers(iconFile,nGRU,indx_meta,err,message) ixHRUfile_min=huge(1) ixHRUfile_max=0 - ! find the min and max hru indices in the state file (why? this should be known from the ) + ! find the min and max hru indices in the state file do iGRU = 1,nGRU do iHRU = 1,gru_struc(iGRU)%hruCount if(gru_struc(iGRU)%hruInfo(iHRU)%hru_nc < ixHRUfile_min) ixHRUfile_min = gru_struc(iGRU)%hruInfo(iHRU)%hru_nc @@ -108,22 +108,18 @@ subroutine read_icond_nlayers(iconFile,nGRU,indx_meta,err,message) end do end do - ! tmp-note: nGRU = size(gru_struc) ... has number of grus in requested gru subset run - ! iHRU_global contains the correct subset index + ! loop over grus in current run to update snow/soil layer information do iGRU = 1,nGRU do iHRU = 1,gru_struc(iGRU)%hruCount iHRU_global = gru_struc(iGRU)%hruInfo(iHRU)%hru_nc - iHRU_local = (iHRU_global - ixHRUfile_min) + 1 - ! single HRU - if(restartFileType==singleHRU)then ! restartFileType hardwired above (to multiHRU) + ! single HRU (Note: 'restartFileType' is hardwired above to multiHRU) + if(restartFileType==singleHRU) then gru_struc(iGRU)%hruInfo(iHRU)%nSnow = snowData(1) gru_struc(iGRU)%hruInfo(iHRU)%nSoil = soilData(1) ! multi HRU else - !gru_struc(iGRU)%hruInfo(iHRU)%nSnow = snowData(iHRU_local) - !gru_struc(iGRU)%hruInfo(iHRU)%nSoil = soilData(iHRU_local) gru_struc(iGRU)%hruInfo(iHRU)%nSnow = snowData(iHRU_global) gru_struc(iGRU)%hruInfo(iHRU)%nSoil = soilData(iHRU_global) endif