Merge pull request #438 from NCAR/develop

arbennett · web-flow · commit 4ee457df3d3c · 2020-11-10T15:44:41.000-08:00
Develop
diff --git 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
@@ -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
@@ -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,14 +117,15 @@ 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
  call read_icond(restartFile,                   & ! intent(in):    name of initial conditions file
                  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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
diff --git a/build/source/netcdf/modelwrite.f90 b/build/source/netcdf/modelwrite.f90
diff --git a/build/source/netcdf/read_icond.f90 b/build/source/netcdf/read_icond.f90
