Continuing PR for SUNDIALS branch

build/source/driver/summa_restart.f90

@@ -66,13 +66,8 @@ subroutine summa_readRestart(summa1_struc, err, message)
  USE globalData,only:elapsedRestart                          ! elapsed time to read model restart files
  ! model decisions
  USE mDecisions_module,only:&                                ! look-up values for the choice of method for the spatial representation of groundwater
-  localColumn,     & ! separate groundwater representation in each local soil column
-  singleBasin        ! single groundwater store over the entire basin
-! look-up values for the numerical method
-USE mDecisions_module,only:&
-  homegrown,       & ! homegrown backward Euler solution using concepts from numerical recipes
-  kinsol,          & ! SUNDIALS backward Euler solution using Kinsol
-  ida                ! SUNDIALS solution using IDA
+   localColumn,    & ! separate groundwater representation in each local soil column
+   singleBasin       ! single groundwater store over the entire basin
  ! look-up values for the choice of variable in energy equations (BE residual or IDA state variable)
  USE mDecisions_module,only:&
    closedForm,     & ! use temperature with closed form heat capacity
@@ -102,7 +97,6 @@ subroutine summa_readRestart(summa1_struc, err, message)
  ! associate to elements in the data structure
  summaVars: associate(& 
   ! model decisions
-  ixNumericalMethod    => model_decisions(iLookDECISIONS%num_method)%iDecision   ,& !choice of numerical solver
   ixNrgConserv         => model_decisions(iLookDECISIONS%nrgConserv)%iDecision   ,& !choice of variable in either energy backward Euler residual or IDA state variable
   spatial_gw           => model_decisions(iLookDECISIONS%spatial_gw)%iDecision   ,& !choice of method for the spatial representation of groundwater
   aquiferIni           => model_decisions(iLookDECISIONS%aquiferIni)%iDecision   ,& !choice of full or empty aquifer at start
@@ -148,7 +142,7 @@ subroutine summa_readRestart(summa1_struc, err, message)
                  progStruct,                    & ! intent(inout): model prognostic variables
                  bvarStruct,                    & ! intent(inout): model basin (GRU) variables
                  indxStruct,                    & ! intent(inout): model indices
-                 no_icond_enth,                 & ! intent(in):    flag that enthalpy not in initial conditions
+                 no_icond_enth,                 & ! intent(out):   flag that enthalpy not in initial conditions
                  err,cmessage)                    ! intent(out):   error control
  if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
 

build/source/dshare/get_ixname.f90

@@ -596,7 +596,7 @@ function get_ixDiag(varName)
   case('scalarVolLatHt_fus'             ); get_ixDiag = iLookDIAG%scalarVolLatHt_fus               ! volumetric latent heat of fusion     (J m-3)
   ! timing information
   case('numFluxCalls'                   ); get_ixDiag = iLookDIAG%numFluxCalls                     ! number of flux calls (-)
-  case('wallClockTime'                  ); get_ixDiag = iLookDIAG%wallClockTime                    ! wall clock time (s)
+  case('wallClockTime'                  ); get_ixDiag = iLookDIAG%wallClockTime                    ! wall clock time for physics routines (s)
   case('meanStepSize'                   ); get_ixDiag = iLookDIAG%meanStepSize                     ! mean time step size (s) over data window
   ! balances
   case('balanceCasNrg'                  ); get_ixDiag = iLookDIAG%balanceCasNrg                    ! balance of energy in the canopy air space (W m-3)
@@ -840,7 +840,7 @@ function get_ixDeriv(varName)
   case('dAquiferTrans_dTGround'         ); get_ixDeriv = iLookDERIV%dAquiferTrans_dTGround         ! derivative in the aquifer transpiration flux w.r.t. ground temperature
   case('dAquiferTrans_dCanWat'          ); get_ixDeriv = iLookDERIV%dAquiferTrans_dCanWat          ! derivative in the aquifer transpiration flux w.r.t. canopy total water
  ! derivative in liquid water fluxes for the soil and snow domain w.r.t temperature
-  case('dFracLiqSnow_dTk'               ); get_ixDeriv = iLookDERIV%dFracLiqSnow_dTk               ! derivative in fraction of liquid snow w.r.t. temperature
+  case('dFracLiqWat_dTk'                ); get_ixDeriv = iLookDERIV%dFracLiqWat_dTk                ! derivative in fraction of liquid water w.r.t. temperature
   case('mLayerdTheta_dTk'               ); get_ixDeriv = iLookDERIV%mLayerdTheta_dTk               ! derivative of volumetric liquid water content w.r.t. temperature (K-1)
   case('mLayerd2Theta_dTk2'             ); get_ixDeriv = iLookDERIV%mLayerd2Theta_dTk2             ! second derivative of volumetric liquid water content w.r.t. temperature
   ! derivatives in time

build/source/dshare/popMetadat.f90

@@ -449,7 +449,7 @@ subroutine popMetadat(err,message)
   diag_meta(iLookDIAG%scalarVolLatHt_fus)              = var_info('scalarVolLatHt_fus'             , 'volumetric latent heat of fusion'                                 , 'J m-3'           , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
   ! timing information
   diag_meta(iLookDIAG%numFluxCalls)                    = var_info('numFluxCalls'                   , 'number of flux calls'                                             , '-'               , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
-  diag_meta(iLookDIAG%wallClockTime)                   = var_info('wallClockTime'                  , 'wall clock time'                                                  , 's'               , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
+  diag_meta(iLookDIAG%wallClockTime)                   = var_info('wallClockTime'                  , 'wall clock time for physics routines'                             , 's'               , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
   diag_meta(iLookDIAG%meanStepSize)                    = var_info('meanStepSize'                   , 'mean time step size over data window'                             , 's'               , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
   ! balances
   diag_meta(iLookDIAG%balanceCasNrg)                   = var_info('balanceCasNrg'                  , 'balance of energy in the canopy air space on data window'         , 'W m-3'           , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
@@ -664,8 +664,8 @@ subroutine popMetadat(err,message)
   deriv_meta(iLookDERIV%dAquiferTrans_dTCanopy)        = var_info('dAquiferTrans_dTCanopy'       , 'derivative in the aquifer transpiration flux w.r.t. canopy temperature',    'm s-1 K-1'  , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
   deriv_meta(iLookDERIV%dAquiferTrans_dTGround)        = var_info('dAquiferTrans_dTGround'       , 'derivative in the aquifer transpiration flux w.r.t. ground temperature',    'm s-1 K-1'  , get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
   deriv_meta(iLookDERIV%dAquiferTrans_dCanWat)         = var_info('dAquiferTrans_dCanWat'        , 'derivative in the aquifer transpiration flux w.r.t. canopy total water',   'm-1 s-1 kg-1', get_ixVarType('scalarv'), iMissVec, iMissVec, .false.)
-  ! derivative in liquid water fluxes for the soil and snow domain w.rt temperatur
-  deriv_meta(iLookDERIV%dFracLiqSnow_dTk)              = var_info('dFracLiqSnow_dTk'             , 'derivative in fraction of liquid snow w.r.t. temperature'             , 'K-1'            , get_ixVarType('midToto'), iMissVec, iMissVec, .false.)
+  ! derivative in liquid water fluxes for the soil and snow domain w.rt temperature
+  deriv_meta(iLookDERIV%dFracLiqWat_dTk)               = var_info('dFracLiqWat_dTk'              , 'derivative in fraction of liquid water w.r.t. temperature'            , 'K-1'            , get_ixVarType('midToto'), iMissVec, iMissVec, .false.)
   deriv_meta(iLookDERIV%mLayerdTheta_dTk)              = var_info('mLayerdTheta_dTk'             , 'derivative of volumetric liquid water content w.r.t. temperature'     , 'K-1'            , get_ixVarType('midToto'), iMissVec, iMissVec, .false.)
   deriv_meta(iLookDERIV%mLayerd2Theta_dTk2)            = var_info('mLayerd2Theta_dTk2'           , 'second derivative of volumetric liquid water content w.r.t. temperature','K-2'           , get_ixVarType('midToto'), iMissVec, iMissVec, .false.)
   ! derivatives in time
@@ -970,12 +970,25 @@ subroutine read_output_file(err,message)
           err=20; return
         endif
 
-      ! temporally constant variables use timestep-level output (no aggregation)
-      case default
-        freqName = trim(lineWords(freqIndex))
-        write(*,*)'WARNING: temporally constant variable '//trim(varName)//': outputting variable in timestep file and will be the same value throughout file'
-        iFreq    = iLookFREQ%timestep
+      ! time and temporally constant variables always outputted at timestep level (no aggregation)
+      case('bpar','attr','type','mpar','time')
+        if(nWords<freqIndex) then
+          freqName = 'empty'
+        else
+          freqName = trim(lineWords(freqIndex))
+        endif
+        if(trim(structName)=='time') then
+          if (freqName/='timestep'.or. freqName/='1') then
+            write(*,*)'WARNING: time variable '//trim(varName)//': outputting variable at timestep level since it cannot be aggregated [entered "'//trim(freqName)//'"]'
+          endif
+        else
+          write(*,*)'WARNING: temporally constant variable '//trim(varName)//': outputting variable in timestep file with no time dimension'
+        endif
+        iFreq = iLookFREQ%timestep
         freqName = 'timestep'
+
+      ! error control
+      case default;  err=20;message=trim(message)//'unable to identify lookup structure';return
     end select
 
     ! --- identify the desired statistic in the metadata structure  -----------

build/source/dshare/var_lookup.f90

@@ -477,7 +477,7 @@ MODULE var_lookup
   integer(i4b)    :: scalarVolLatHt_fus              = integerMissing ! volumetric latent heat of fusion     (J m-3)
   ! number of function evaluations
   integer(i4b)    :: numFluxCalls                    = integerMissing ! number of flux calls (-)
-  integer(i4b)    :: wallClockTime                   = integerMissing ! wall clock time (s)
+  integer(i4b)    :: wallClockTime                   = integerMissing ! wall clock time for physics routines(s)
   integer(i4b)    :: meanStepSize                    = integerMissing ! mean time step size over data window (s)
   ! balances
   integer(i4b)    :: balanceCasNrg                   = integerMissing ! balance of energy in the canopy air space (W m-3)
@@ -698,7 +698,7 @@ MODULE var_lookup
   integer(i4b)    :: dAquiferTrans_dTGround          = integerMissing ! derivative in the aquifer transpiration flux w.r.t. ground temperature
   integer(i4b)    :: dAquiferTrans_dCanWat           = integerMissing ! derivative in the aquifer transpiration flux w.r.t. canopy total water
  ! derivative in liquid water fluxes for the soil and snow domain w.r.t temperature
-  integer(i4b)    :: dFracLiqSnow_dTk                = integerMissing ! derivative in fraction of liquid snow w.r.t. temperature
+  integer(i4b)    :: dFracLiqWat_dTk                 = integerMissing ! derivative in fraction of liquid water w.r.t. temperature
   integer(i4b)    :: mLayerdTheta_dTk                = integerMissing ! derivative of volumetric liquid water content w.r.t. temperature (K-1)
   integer(i4b)    :: mLayerd2Theta_dTk2              = integerMissing ! second derivative of volumetric liquid water content w.r.t. temperature
  ! derivatives in time

build/source/engine/computHeatCap.f90

@@ -339,7 +339,6 @@ end subroutine computHeatCapAnalytic
 !   NOTE: computing on whole vector, could just compute on state subset
 ! **********************************************************************************************************
 subroutine computCm(&
-                      be_solver,               & ! intent(in):  flag for BE solver, need to include latent heat part in Jacobian Cm if true
                       ! input: state variables
                       canopyDepth,             & ! intent(in):  depth of the vegetation canopy (m)
                       scalarCanopyTemp,        & ! intent(in):  value of canopy temperature (K)
@@ -350,9 +349,7 @@ subroutine computCm(&
                       indx_data,               & ! intent(in):  model layer indices
                       ! output
                       scalarCanopyCm,          & ! intent(inout): Cm for vegetation (J kg K-1)
-                      scalarCanopyCm_noLH,     & ! intent(inout): Cm without latent heat part for vegetation (J kg K-1)
-                      mLayerCm,                & ! intent(inout): Cm for soil and snow (J kg K-1)
-                      mLayerCm_noLH,           & ! intent(inout): Cm without latent heat part for soil and snow (J kg K-1)
+                      mLayerCm,                & ! intent(inout): Cm for snow and soil (J kg K-1)
                       dCm_dPsi0,               & ! intent(inout): derivative in Cm w.r.t. matric potential (J kg)
                       dCm_dTk,                 & ! intent(inout): derivative in Cm w.r.t. temperature (J kg K-2)
                       dCm_dTkCanopy,           & ! intent(inout): derivative in Cm w.r.t. temperature (J kg K-2)
@@ -365,7 +362,6 @@ subroutine computCm(&
   USE soil_utils_module,only:crit_soilT     ! compute critical temperature below which ice exists (soil)
   ! --------------------------------------------------------------------------------------------------------------------------------------
   ! input: state variables
-  logical(i4b),intent(in)              :: be_solver              ! flag for BE solver, need to include latent heat part in Jacobian Cm if true
   real(rkind),intent(in)               :: canopyDepth            ! depth of the vegetation canopy (m)
   real(rkind),intent(in)               :: scalarCanopyTemp       ! value of canopy temperature (K)
   real(rkind),intent(in)               :: mLayerTemp(:)          ! vector of temperature (K)
@@ -375,9 +371,7 @@ subroutine computCm(&
   type(var_ilength),intent(in)         :: indx_data              ! model layer indices
   ! output: Cm and derivatives
   real(rkind),intent(inout)            :: scalarCanopyCm         ! Cm for vegetation (J kg K-1) use for LHS
-  real(rkind),intent(inout)            :: scalarCanopyCm_noLH    ! Cm without latent heat part for vegetation (J kg K-1) use for RHS
-  real(rkind),intent(inout)            :: mLayerCm(:)            ! Cm for soil and snow (J kg K-1)
-  real(rkind),intent(inout)            :: mLayerCm_noLH(:)       ! Cm without latent heat part for soil and snow (J kg K-1)
+  real(rkind),intent(inout)            :: mLayerCm(:)            ! Cm for snow and soil (J kg K-1)
   real(rkind),intent(inout)            :: dCm_dPsi0(:)           ! derivative in Cm w.r.t. matric potential (J kg)
   real(rkind),intent(inout)            :: dCm_dTk(:)             ! derivative in Cm w.r.t. temperature (J kg K-2)
   real(rkind),intent(inout)            :: dCm_dTkCanopy          ! derivative in Cm w.r.t. temperature (J kg K-2)
@@ -448,52 +442,42 @@ subroutine computCm(&
             ! Note that scalarCanopyCm/iden_water is computed
             diffT = scalarCanopyTemp - Tfreeze
             if(diffT>=0._rkind)then
-              scalarCanopyCm_noLH =  Cp_water * diffT
-              scalarCanopyCm      =  scalarCanopyCm_noLH
+              scalarCanopyCm =  Cp_water * diffT
               ! derivatives
               dCm_dTkCanopy  = Cp_water
             else
               integral = (1._rkind/snowfrz_scale) * atan(snowfrz_scale * diffT)
               fLiq = fracLiquid(scalarCanopyTemp,snowfrz_scale)
-              scalarCanopyCm_noLH = Cp_water * integral + Cp_ice * (diffT - integral) 
-              scalarCanopyCm      = scalarCanopyCm_noLH
-              if (be_solver) scalarCanopyCm = scalarCanopyCm_noLH - LH_fus * (1._rkind - fLiq) / canopyDepth ! LH_fus is term is not the Jacobian for the RHS
+              scalarCanopyCm = Cp_water * integral + Cp_ice * (diffT - integral) 
               ! derivatives
               dfLiq_dT = dFracLiq_dTk(scalarCanopyTemp,snowfrz_scale)
               dCm_dTkCanopy = Cp_water * fLiq + Cp_ice * (1._rkind - fLiq)
-              if (be_solver) dCm_dTkCanopy = dCm_dTkCanopy + LH_fus * dfLiq_dT / canopyDepth ! LH_fus is term is not in the Jacobian for the LHS
             end if
 
           case(iname_snow)
             diffT = mLayerTemp(iLayer) - Tfreeze
             fLiq = fracLiquid(mLayerTemp(iLayer),snowfrz_scale)
             integral = (1._rkind/snowfrz_scale) * atan(snowfrz_scale * diffT)
-            mLayerCm_noLH(iLayer) = (iden_water * Cp_ice - iden_air * Cp_air * iden_water/iden_ice) * ( diffT - integral ) &
+            mLayerCm(iLayer) = (iden_water * Cp_ice - iden_air * Cp_air * iden_water/iden_ice) * ( diffT - integral ) &
                                    + (iden_water * Cp_water - iden_air * Cp_air) * integral
-            mLayerCm(iLayer)      = mLayerCm_noLH(iLayer)
-            if (be_solver) mLayerCm(iLayer) = mLayerCm_noLH(iLayer) - LH_fus * iden_water * (1._rkind - fLiq) ! LH_fus is term is already in the Jacobian for the LHS
             ! derivatives
             dfLiq_dT = dFracLiq_dTk(mLayerTemp(iLayer),snowfrz_scale)
             dCm_dTk(iLayer) = (iden_water * Cp_ice - iden_air * Cp_air * iden_water/iden_ice) * ( 1._rkind -fLiq ) &
                              + (iden_water * Cp_water - iden_air * Cp_air) * fLiq
-            if (be_solver) dCm_dTk(iLayer) = dCm_dTk(iLayer) + LH_fus * iden_water * dfLiq_dT ! LH_fus is term is not in the Jacobian for the LHS
 
           case(iname_soil)
             diffT = mLayerTemp(iLayer) - Tfreeze
             Tcrit = crit_soilT( mLayerMatricHead(ixControlIndex) )
             diff0 = Tcrit - Tfreeze
             if( mLayerTemp(iLayer)>=Tcrit)then
-              mLayerCm_noLH(iLayer) = (-iden_air * Cp_air + iden_water * Cp_water) * diffT
-              mLayerCm(iLayer)      = mLayerCm_noLH(iLayer)
+              mLayerCm(iLayer) = (-iden_air * Cp_air + iden_water * Cp_water) * diffT
               ! derivatives
               dCm_dTk(iLayer) = -iden_air * Cp_air + iden_water * Cp_water
               dCm_dPsi0(ixControlIndex) = 0._rkind
             else        
-              mLayerCm_noLH(iLayer) = -iden_air * Cp_air * diffT + iden_ice * Cp_ice * (mLayerTemp(iLayer)-Tcrit) &
+              mLayerCm(iLayer) = -iden_air * Cp_air * diffT + iden_ice * Cp_ice * (mLayerTemp(iLayer)-Tcrit) &
                                      + iden_water * Cp_water * diff0
-              mLayerCm(iLayer)      = mLayerCm_noLH(iLayer)
-              if (be_solver) mLayerCm(iLayer) = mLayerCm_noLH(iLayer) - LH_fus * iden_water ! LH_fus is term is already in the Jacobian for the LHS
-              ! derivatives, note that does not matter if be_solver is true or false
+              ! derivatives
               dTcrit_dPsi0 = merge(gravity*Tfreeze/LH_fus,0._rkind,mLayerMatricHead(ixControlIndex)<=0._rkind)
               dCm_dTk(iLayer) = -iden_air * Cp_air + iden_ice * Cp_ice
               dCm_dPsi0(ixControlIndex) = (-iden_ice * Cp_ice + iden_water * Cp_water) * dTcrit_dPsi0