Merge pull request #560 from ashleymedin/master

Allocation error fixes

build/source/engine/eval8summa.f90

@@ -483,7 +483,7 @@ subroutine eval8summa(&
  endif
 
  ! snow+soil domain: get the correct water states (total water, or liquid water, depending on the state type)
- mLayerVolFracHydTrial = merge(mLayerVolFracWatTrial, mLayerVolFracLiqTrial, (ixHydType==iname_watLayer .or. ixHydType==iname_matLayer) )
+ mLayerVolFracHydTrial = merge(mLayerVolFracWatTrial, mLayerVolFracLiqTrial, (ixHydType(1:nLayers)==iname_watLayer .or. ixHydType(1:nLayers)==iname_matLayer) )
 
  ! compute the residual vector
  call computResid(&

build/source/engine/groundwatr.f90

@@ -438,7 +438,7 @@ subroutine computeBaseflow(&
  totalColumnOutflow = sum(mLayerColumnOutflow(1:nSoil))/HRUarea
 
  ! compute the available storage (m)
- availStorage = sum(mLayerDepth(1:nSoil)*(theta_sat - (mLayerVolFracLiq(1:nSoil)+mLayerVolFracIce(1:nSoil))) )
+ availStorage = sum(mLayerDepth(1:nSoil)*(theta_sat(1:nSoil) - (mLayerVolFracLiq(1:nSoil)+mLayerVolFracIce(1:nSoil))))
 
  ! compute the smoothing function (-)
  if(availStorage < xMinEval)then

build/source/engine/indexState.f90

@@ -381,7 +381,7 @@ subroutine indexSplit(stateSubsetMask,             & ! intent(in)    : logical v
  ! -------------------------------------------
 
  ! get different masks
- volFracWat_mask = (ixHydType==iname_watLayer .or. ixHydType==iname_liqLayer)
+ volFracWat_mask = (ixHydType(      1:nLayers)==iname_watLayer .or. ixHydType(      1:nLayers)==iname_liqLayer)
  matricHead_mask = (ixHydType(nSnow+1:nLayers)==iname_matLayer .or. ixHydType(nSnow+1:nLayers)==iname_lmpLayer)
 
  ! get state subsets for desired variables

build/source/engine/layerDivide.f90

@@ -330,8 +330,8 @@ subroutine layerDivide(&
   layerType(nSnow+2:nLayers+1) = iname_soil
 
   ! identify the number of snow and soil layers, and check all is a-OK
-  nSnow   = count(layerType==iname_snow)
-  nSoil   = count(layerType==iname_soil)
+  nSnow   = count(layerType(1:nLayers+1)==iname_snow)
+  nSoil   = count(layerType(1:nLayers+1)==iname_soil)
   nLayers = nSnow + nSoil
 
   ! re-set coordinate variables

build/source/engine/opSplittin.f90

@@ -1110,6 +1110,9 @@ subroutine stateFilter(ixCoupling,ixSolution,ixStateThenDomain,iStateTypeSplit,i
  ! identify scalar solutions
  if(ixSolution==scalar)then
 
+  allocate(ixSubset(1_i4b))
+  ixSubset = 0._rkind
+
   ! get the subset of indices
   call indxSubset(ixSubset, ixAllState, stateMask, err, cmessage)
   if(err/=0)then; message=trim(message)//trim(cmessage); return; endif

build/source/engine/paramCheck.f90

@@ -116,8 +116,8 @@ subroutine paramCheck(mpar_data,err,message)
  heightCanopyTop        => mpar_data%var(iLookPARAM%heightCanopyTop)%dat(1),   & ! intent(in): [dp] height at the top of the vegetation canopy (m)
  heightCanopyBottom     => mpar_data%var(iLookPARAM%heightCanopyBottom)%dat(1),& ! intent(in): [dp] height at the bottom of the vegetation canopy (m)
  ! transpiration
- critSoilWilting        => mpar_data%var(iLookPARAM%critSoilWilting)%dat,      & ! intent(in): [dp] critical vol. liq. water content when plants are wilting (-)
- critSoilTranspire      => mpar_data%var(iLookPARAM%critSoilTranspire)%dat,    & ! intent(in): [dp] critical vol. liq. water content when transpiration is limited (-)
+ critSoilWilting        => mpar_data%var(iLookPARAM%critSoilWilting)%dat(1),   & ! intent(in): [dp] critical vol. liq. water content when plants are wilting (-)
+ critSoilTranspire      => mpar_data%var(iLookPARAM%critSoilTranspire)%dat(1), & ! intent(in): [dp] critical vol. liq. water content when transpiration is limited (-)
  ! soil properties
  fieldCapacity          => mpar_data%var(iLookPARAM%fieldCapacity)%dat(1),     & ! intent(in): [dp]    field capacity (-)
  theta_sat              => mpar_data%var(iLookPARAM%theta_sat)%dat,            & ! intent(in): [dp(:)] soil porosity (-)
@@ -161,7 +161,7 @@ subroutine paramCheck(mpar_data,err,message)
  end if
 
  ! check transpiration
- if( any(critSoilTranspire < critSoilWilting) )then
+ if( critSoilTranspire < critSoilWilting )then
   write(message,'(a,i0,a)') trim(message)//'critical point for transpiration is less than the wilting point'
   err=20; return
  endif

build/source/engine/summaSolve.f90

@@ -375,6 +375,7 @@ subroutine lineSearchRefinement(doLineSearch,stateVecTrial,newtStepScaled,aJacSc
   ! --------------------------------------------------------------------------------------------------------
   ! initialize error control
   err=0; message='lineSearchRefinement/'
+  converged =.false.
 
   ! check the need to compute the line search
   if(doLineSearch)then
@@ -531,6 +532,7 @@ subroutine trustRegionRefinement(doTrustRefinement,stateVecTrial,newtStepScaled,
 
   ! --------------------------------------------------------------------------------------------------------
   err=0; message='trustRegionRefinement/'
+  converged =.false.
 
   ! check the need to refine the step
   if(doTrustRefinement)then
@@ -603,6 +605,7 @@ subroutine safeRootfinder(stateVecTrial,rVecscaled,newtStepScaled,stateVecNew,fl
   !real(rkind)                      :: xIncrement(nState)       ! trial increment (not used)
   ! --------------------------------------------------------------------------------------------------------
   err=0; message='safeRootfinder/'
+  converged = .false.
 
   ! check scalar
   if(size(stateVecTrial)/=1 .or. size(rVecScaled)/=1 .or. size(newtStepScaled)/=1)then

build/source/engine/var_derive.f90

@@ -462,7 +462,7 @@ subroutine fracFuture(bpar_data,bvar_data,err,message)
    end do ! (looping through future time steps)
 
    ! check that we have enough bins
-   sumFrac  = sum(fractionFuture)
+   sumFrac  = sum(fractionFuture(1:nTDH))
    if(abs(1._rkind - sumFrac) > tolerFrac)then
     write(*,*) 'WARNING: The fraction of basin runoff histogram being accounted for by time delay vector is ', sumFrac
     write(*,*) 'This is less than allowed by tolerFrac = ', tolerFrac
@@ -473,7 +473,7 @@ subroutine fracFuture(bpar_data,bvar_data,err,message)
     write(*,*) '       -- note that nTimeDelay defines the number of time steps in the time delay histogram'
    end if
    ! ensure the fraction sums to one
-   fractionFuture = fractionFuture/sumFrac
+   fractionFuture(1:nTDH) = fractionFuture(1:nTDH)/sumFrac
 
   ! ** error checking
   case default; err=20; message=trim(message)//'cannot find option for sub-grid routing'; return