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..61b80816e 100755 --- a/build/source/driver/summa_restart.f90 +++ b/build/source/driver/summa_restart.f90 @@ -101,6 +101,7 @@ 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/' @@ -116,7 +117,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 +125,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 +145,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 +180,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 +227,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..32ef3d466 100755 --- a/build/source/driver/summa_writeOutput.f90 +++ b/build/source/driver/summa_writeOutput.f90 @@ -287,7 +287,7 @@ 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/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/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/dshare/globalData.f90 b/build/source/dshare/globalData.f90 index 5fb4301ac..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 @@ -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 + 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 e7e41ed5d..27e73300a 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 -integer(i4b),parameter :: nBand=2 ! number of spectral bands -integer(i4b),parameter :: 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,7 +544,8 @@ 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 @@ -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 ba61a20f1..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 @@ -272,8 +273,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 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/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/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/engine/var_derive.f90 b/build/source/engine/var_derive.f90 index cdf93f990..8227b0407 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,10 +415,10 @@ 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 + ! initialize runoffFuture (will be overwritten by initial conditions file values if present) runoffFuture(1:nTDH) = 0._dp ! select option for sub-grid routing @@ -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 globalData.f90' err=20; return end if ! ensure the fraction sums to one diff --git a/build/source/netcdf/def_output.f90 b/build/source/netcdf/def_output.f90 index 792196735..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,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)%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 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 diff --git a/build/source/netcdf/modelwrite.f90 b/build/source/netcdf/modelwrite.f90 index dc417a9f4..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 @@ -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,24 @@ 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 globalData,only:nTimeDelay ! number of timesteps in the time delay 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 @@ -488,8 +495,11 @@ subroutine writeRestart(filename, & ! intent(in): name of restart file integer(i4b) :: nLayers ! number of total layers 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,37 @@ 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+1)) ! include 1 additional basin variable in ID array (possibly more later) ! maximum number of soil layers maxSoil = gru_struc(1)%hruInfo(1)%nSoil ! maximum number of snow layers maxSnow = maxSnowLayers - + ! 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) ,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 + 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 +592,11 @@ 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) ! 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 +664,18 @@ 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,nTimeDelay/)) + + 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 old mode 100755 new mode 100644 index 6559af4a7..c5bdd929e --- 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 @@ -41,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 @@ -78,29 +80,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 sotrage 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 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,20 +108,20 @@ subroutine read_icond_nlayers(iconFile,nGRU,indx_meta,err,message) end do end do + ! 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 + ! 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 end do @@ -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,19 +153,24 @@ 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: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 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 USE updatState_module,only:updateSoil ! update soil states + implicit none ! -------------------------------------------------------------------------------------------------------- @@ -174,6 +180,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,7 +188,9 @@ 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(1) :: ndx ! intermediate array of loop indices integer(i4b) :: iGRU ! loop index integer(i4b) :: iHRU ! loop index integer(i4b) :: dimID ! varible dimension ids @@ -194,13 +203,16 @@ 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),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 ! -------------------------------------------------------------------------------------------------------- @@ -256,21 +268,22 @@ 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) - 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 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 @@ -279,10 +292,9 @@ 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 - ! 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 @@ -319,9 +331,9 @@ 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 + if(err/=0)then; message=trim(message)//'problem deallocating HRU variable data'; return; endif - end do ! iVar + end do ! end looping through prognostic variables (iVar) ! -------------------------------------------------------------------------------------------------------- ! (2) set number of layers @@ -342,7 +354,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 @@ -373,10 +385,82 @@ 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 + ! -------------------------------------------------------------------------------------------------------- + ! (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); + 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 + 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 (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 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+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 + + 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 + end subroutine read_icond end module read_icond_module