From e962c7b031af3aeeecb1db8a80ab2f25614d1a1a Mon Sep 17 00:00:00 2001 From: Kyle <kyle.c.klenk@gmail.com> Date: Fri, 9 Sep 2022 19:07:56 +0000 Subject: [PATCH] changed jacob and resid --- build/source/engine/sundials/computJacDAE.f90 | 736 ----------- .../engine/sundials/computJacobSundials.f90 | 1120 +++++++++++++++++ .../source/engine/sundials/computResidDAE.f90 | 346 ----- .../engine/sundials/computResidSundials.f90 | 347 +++++ utils/laugh_tests/colbeck1976/actors_out.txt | 349 +++++ 5 files changed, 1816 insertions(+), 1082 deletions(-) delete mode 100644 build/source/engine/sundials/computJacDAE.f90 create mode 100644 build/source/engine/sundials/computJacobSundials.f90 delete mode 100644 build/source/engine/sundials/computResidDAE.f90 create mode 100644 build/source/engine/sundials/computResidSundials.f90 create mode 100644 utils/laugh_tests/colbeck1976/actors_out.txt diff --git a/build/source/engine/sundials/computJacDAE.f90 b/build/source/engine/sundials/computJacDAE.f90 deleted file mode 100644 index f1f87d2..0000000 --- a/build/source/engine/sundials/computJacDAE.f90 +++ /dev/null @@ -1,736 +0,0 @@ -! SUMMA - Structure for Unifying Multiple Modeling Alternatives -! Copyright (C) 2014-2015 NCAR/RAL -! -! This file is part of SUMMA -! -! For more information see: http://www.ral.ucar.edu/projects/summa -! -! This program is free software: you can redistribute it and/or modify -! it under the terms of the GNU General Public License as published by -! the Free Software Foundation, either version 3 of the License, or -! (at your option) any later version. -! -! This program is distributed in the hope that it will be useful, -! but WITHOUT ANY WARRANTY; without even the implied warranty of -! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -! GNU General Public License for more details. -! -! You should have received a copy of the GNU General Public License -! along with this program. If not, see <http://www.gnu.org/licenses/>. - -module computJacDAE_module - -! data types -USE nrtype - -! derived types to define the data structures -USE data_types,only:& - var_ilength, & ! data vector with variable length dimension (i4b) - var_dlength ! data vector with variable length dimension (rkind) - -! named variables for structure elements -USE var_lookup,only:iLookPROG ! named variables for structure elements -USE var_lookup,only:iLookDIAG ! named variables for structure elements -USE var_lookup,only:iLookINDEX ! named variables for structure elements -USE var_lookup,only:iLookDERIV ! named variables for structure elements - -! look-up values for the form of Richards' equation -USE mDecisions_module,only: & - moisture, & ! moisture-based form of Richards' equation - mixdform ! mixed form of Richards' equation - -! access the global print flag -USE globalData,only:globalPrintFlag - -! access missing values -USE globalData,only:integerMissing ! missing integer -USE globalData,only:realMissing ! missing real number - -! domain types -USE globalData,only:iname_veg ! named variables for vegetation -USE globalData,only:iname_snow ! named variables for snow -USE globalData,only:iname_soil ! named variables for soil - -! named variables to describe the state variable type -USE globalData,only:iname_nrgCanair ! named variable defining the energy of the canopy air space -USE globalData,only:iname_nrgCanopy ! named variable defining the energy of the vegetation canopy -USE globalData,only:iname_watCanopy ! named variable defining the mass of water on the vegetation canopy -USE globalData,only:iname_nrgLayer ! named variable defining the energy state variable for snow+soil layers -USE globalData,only:iname_watLayer ! named variable defining the total water state variable for snow+soil layers -USE globalData,only:iname_liqLayer ! named variable defining the liquid water state variable for snow+soil layers -USE globalData,only:iname_matLayer ! named variable defining the matric head state variable for soil layers -USE globalData,only:iname_lmpLayer ! named variable defining the liquid matric potential state variable for soil layers - -! access named variables to describe the form and structure of the matrices used in the numerical solver -USE globalData,only: ku ! number of super-diagonal bands -USE globalData,only: kl ! number of sub-diagonal bands -USE globalData,only: ixDiag ! index for the diagonal band -USE globalData,only: nBands ! length of the leading dimension of the band diagonal matrix -USE globalData,only: ixFullMatrix ! named variable for the full Jacobian matrix -USE globalData,only: ixBandMatrix ! named variable for the band diagonal matrix -USE globalData,only: iJac1 ! first layer of the Jacobian to print -USE globalData,only: iJac2 ! last layer of the Jacobian to print - -! constants -USE multiconst,only:& - LH_fus, & ! latent heat of fusion (J kg-1) - iden_water, & ! intrinsic density of liquid water (kg m-3) - ! specific heat - Cp_ice, & ! specific heat of ice (J kg-1 K-1) - Cp_water ! specific heat of liquid water (J kg-1 K-1) - -implicit none -! define constants -real(rkind),parameter :: verySmall=tiny(1.0_rkind) ! a very small number -integer(i4b),parameter :: ixBandOffset=kl+ku+1 ! offset in the band Jacobian matrix - -private -public::computJacDAE -contains - - ! ********************************************************************************************************** - ! public subroutine computJacDAE: compute the Jacobian matrix - ! ********************************************************************************************************** - subroutine computJacDAE(& - ! input: model control - cj, & ! intent(in): this scalar changes whenever the step size or method order changes - dt, & ! intent(in): length of the time step (seconds) - nSnow, & ! intent(in): number of snow layers - nSoil, & ! intent(in): number of soil layers - nLayers, & ! intent(in): total number of layers - computeVegFlux, & ! intent(in): flag to indicate if we need to compute fluxes over vegetation - computeBaseflow, & ! intent(in): flag to indicate if we need to compute baseflow - ixMatrix, & ! intent(in): form of the Jacobian matrix - specificStorage, & ! intent(in): specific storage coefficient (m-1) - theta_sat, & ! intent(in): soil porosity (-) - ixRichards, & ! intent(in): choice of option for Richards' equation - ! input: data structures - indx_data, & ! intent(in): index data - prog_data, & ! intent(in): model prognostic variables for a local HRU - diag_data, & ! intent(in): model diagnostic variables for a local HRU - deriv_data, & ! intent(in): derivatives in model fluxes w.r.t. relevant state variables - dBaseflow_dMatric, & ! intent(in): derivative in baseflow w.r.t. matric head (s-1) - ! input: state variables - mLayerTemp, & ! intent(in): vector of layer temperature (K) - mLayerTempPrime, & ! intent(in) - mLayerMatricHeadPrime, & ! intent(in) - mLayerMatricHeadLiqPrime, & ! intent(in) - mLayerVolFracWatPrime, & ! intent(in) - scalarCanopyTemp, & ! intent(in) - scalarCanopyTempPrime, & ! intent(in) derivative value for temperature of the vegetation canopy (K) - scalarCanopyWatPrime, & ! intent(in) - ! input-output: Jacobian and its diagonal - dMat, & ! intent(inout): diagonal of the Jacobian matrix - aJac, & ! intent(out): Jacobian matrix - ! output: error control - err,message) ! intent(out): error code and error message - ! ----------------------------------------------------------------------------------------------------------------- - implicit none - ! input: model control - real(rkind),intent(in) :: cj - real(rkind),intent(in) :: dt ! length of the time step (seconds) - integer(i4b),intent(in) :: nSnow ! number of snow layers - integer(i4b),intent(in) :: nSoil ! number of soil layers - integer(i4b),intent(in) :: nLayers ! total number of layers in the snow+soil domain - logical(lgt),intent(in) :: computeVegFlux ! flag to indicate if computing fluxes over vegetation - logical(lgt),intent(in) :: computeBaseflow ! flag to indicate if computing baseflow - integer(i4b),intent(in) :: ixMatrix ! form of the Jacobian matrix - real(rkind),intent(in) :: specificStorage ! specific storage coefficient (m-1) - real(rkind),intent(in) :: theta_sat(:) ! soil porosity (-) - integer(i4b),intent(in) :: ixRichards ! choice of option for Richards' equation - ! input: data structures - type(var_ilength),intent(in) :: indx_data ! indices defining model states and layers - type(var_dlength),intent(in) :: prog_data ! prognostic variables for a local HRU - type(var_dlength),intent(in) :: diag_data ! diagnostic variables for a local HRU - type(var_dlength),intent(in) :: deriv_data ! derivatives in model fluxes w.r.t. relevant state variables - real(rkind),intent(in) :: dBaseflow_dMatric(:,:) ! derivative in baseflow w.r.t. matric head (s-1) - ! input: state variables - real(rkind),intent(in) :: mLayerTemp(:) - real(rkind),intent(in) :: mLayerTempPrime(:) - real(rkind),intent(in) :: mLayerMatricHeadPrime(:) - real(rkind),intent(in) :: mLayerMatricHeadLiqPrime(:) - real(rkind),intent(in) :: mLayerVolFracWatPrime(:) - real(rkind),intent(in) :: scalarCanopyTemp - real(rkind),intent(in) :: scalarCanopyTempPrime ! derivative value for temperature of the vegetation canopy (K) - real(rkind),intent(in) :: scalarCanopyWatPrime - - ! input-output: Jacobian and its diagonal - real(rkind),intent(inout) :: dMat(:) ! diagonal of the Jacobian matrix - real(rkind),intent(out) :: aJac(:,:) ! Jacobian matrix - ! output variables - integer(i4b),intent(out) :: err ! error code - character(*),intent(out) :: message ! error message - ! -------------------------------------------------------------- - ! * local variables - ! -------------------------------------------------------------- - ! indices of model state variables - integer(i4b) :: jState ! index of state within the state subset - integer(i4b) :: qState ! index of cross-derivative state variable for baseflow - integer(i4b) :: nrgState ! energy state variable - integer(i4b) :: watState ! hydrology state variable - integer(i4b) :: nState ! number of state variables - ! indices of model layers - integer(i4b) :: iLayer ! index of model layer - integer(i4b) :: jLayer ! index of model layer within the full state vector (hydrology) - integer(i4b) :: pLayer ! indices of soil layers (used for the baseflow derivatives) - ! conversion factors - real(rkind) :: convLiq2tot ! factor to convert liquid water derivative to total water derivative - - ! -------------------------------------------------------------- - ! associate variables from data structures - associate(& - ! indices of model state variables - ixCasNrg => indx_data%var(iLookINDEX%ixCasNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy air space energy state variable - ixVegNrg => indx_data%var(iLookINDEX%ixVegNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy energy state variable - ixVegHyd => indx_data%var(iLookINDEX%ixVegHyd)%dat(1) ,& ! intent(in): [i4b] index of canopy hydrology state variable (mass) - ixTopNrg => indx_data%var(iLookINDEX%ixTopNrg)%dat(1) ,& ! intent(in): [i4b] index of upper-most energy state in the snow+soil subdomain - ixTopHyd => indx_data%var(iLookINDEX%ixTopHyd)%dat(1) ,& ! intent(in): [i4b] index of upper-most hydrology state in the snow+soil subdomain - ixAqWat => indx_data%var(iLookINDEX%ixAqWat)%dat(1) ,& ! intent(in): [i4b] index of water storage in the aquifer - ! vectors of indices for specfic state types within specific sub-domains IN THE FULL STATE VECTOR - ixNrgLayer => indx_data%var(iLookINDEX%ixNrgLayer)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in the snow+soil domain - ixHydLayer => indx_data%var(iLookINDEX%ixHydLayer)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for hydrology states in the snow+soil domain - ! vector of energy indices for the snow and soil domains - ! NOTE: states not in the subset are equal to integerMissing - ixSnowSoilNrg => indx_data%var(iLookINDEX%ixSnowSoilNrg)%dat ,& ! intent(in): [i4b(:)] index in the state subset for energy state variables in the snow+soil domain - ixSnowOnlyNrg => indx_data%var(iLookINDEX%ixSnowOnlyNrg)%dat ,& ! intent(in): [i4b(:)] index in the state subset for energy state variables in the snow domain - ixSoilOnlyNrg => indx_data%var(iLookINDEX%ixSoilOnlyNrg)%dat ,& ! intent(in): [i4b(:)] index in the state subset for energy state variables in the soil domain - ! vector of hydrology indices for the snow and soil domains - ! NOTE: states not in the subset are equal to integerMissing - ixSnowSoilHyd => indx_data%var(iLookINDEX%ixSnowSoilHyd)%dat ,& ! intent(in): [i4b(:)] index in the state subset for hydrology state variables in the snow+soil domain - ixSnowOnlyHyd => indx_data%var(iLookINDEX%ixSnowOnlyHyd)%dat ,& ! intent(in): [i4b(:)] index in the state subset for hydrology state variables in the snow domain - ixSoilOnlyHyd => indx_data%var(iLookINDEX%ixSoilOnlyHyd)%dat ,& ! intent(in): [i4b(:)] index in the state subset for hydrology state variables in the soil domain - ! number of state variables of a specific type - nSnowSoilNrg => indx_data%var(iLookINDEX%nSnowSoilNrg )%dat(1) ,& ! intent(in): [i4b] number of energy state variables in the snow+soil domain - nSnowOnlyNrg => indx_data%var(iLookINDEX%nSnowOnlyNrg )%dat(1) ,& ! intent(in): [i4b] number of energy state variables in the snow domain - nSoilOnlyNrg => indx_data%var(iLookINDEX%nSoilOnlyNrg )%dat(1) ,& ! intent(in): [i4b] number of energy state variables in the soil domain - nSnowSoilHyd => indx_data%var(iLookINDEX%nSnowSoilHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the snow+soil domain - nSnowOnlyHyd => indx_data%var(iLookINDEX%nSnowOnlyHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the snow domain - nSoilOnlyHyd => indx_data%var(iLookINDEX%nSoilOnlyHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the soil domain - ! type and index of model control volume - ixHydType => indx_data%var(iLookINDEX%ixHydType)%dat ,& ! intent(in): [i4b(:)] index of the type of hydrology states in snow+soil domain - ixDomainType => indx_data%var(iLookINDEX%ixDomainType)%dat ,& ! intent(in): [i4b(:)] indices defining the type of the domain (iname_veg, iname_snow, iname_soil) - ixControlVolume => indx_data%var(iLookINDEX%ixControlVolume)%dat ,& ! intent(in): [i4b(:)] index of the control volume for specific model domains - ! mapping between states and model layers - ixMapSubset2Full => indx_data%var(iLookINDEX%ixMapSubset2Full)%dat ,& ! intent(in): [i4b(:)] list of indices in the full state vector that are in the state subset - ixMapFull2Subset => indx_data%var(iLookINDEX%ixMapFull2Subset)%dat ,& ! intent(in): [i4b(:)] list of indices in the state subset in each element of the full state vector - ! derivatives in net vegetation energy fluxes w.r.t. relevant state variables - dCanairNetFlux_dCanairTemp => deriv_data%var(iLookDERIV%dCanairNetFlux_dCanairTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy air space flux w.r.t. canopy air temperature - dCanairNetFlux_dCanopyTemp => deriv_data%var(iLookDERIV%dCanairNetFlux_dCanopyTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy air space flux w.r.t. canopy temperature - dCanairNetFlux_dGroundTemp => deriv_data%var(iLookDERIV%dCanairNetFlux_dGroundTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy air space flux w.r.t. ground temperature - dCanopyNetFlux_dCanairTemp => deriv_data%var(iLookDERIV%dCanopyNetFlux_dCanairTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy flux w.r.t. canopy air temperature - dCanopyNetFlux_dCanopyTemp => deriv_data%var(iLookDERIV%dCanopyNetFlux_dCanopyTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy flux w.r.t. canopy temperature - dCanopyNetFlux_dGroundTemp => deriv_data%var(iLookDERIV%dCanopyNetFlux_dGroundTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy flux w.r.t. ground temperature - dCanopyNetFlux_dCanWat => deriv_data%var(iLookDERIV%dCanopyNetFlux_dCanWat )%dat(1) ,& ! intent(in): [dp] derivative in net canopy fluxes w.r.t. canopy total water content - dGroundNetFlux_dCanairTemp => deriv_data%var(iLookDERIV%dGroundNetFlux_dCanairTemp )%dat(1) ,& ! intent(in): [dp] derivative in net ground flux w.r.t. canopy air temperature - dGroundNetFlux_dCanopyTemp => deriv_data%var(iLookDERIV%dGroundNetFlux_dCanopyTemp )%dat(1) ,& ! intent(in): [dp] derivative in net ground flux w.r.t. canopy temperature - dGroundNetFlux_dCanWat => deriv_data%var(iLookDERIV%dGroundNetFlux_dCanWat )%dat(1) ,& ! intent(in): [dp] derivative in net ground fluxes w.r.t. canopy total water content - ! derivatives in evaporative fluxes w.r.t. relevant state variables - dCanopyEvaporation_dTCanair => deriv_data%var(iLookDERIV%dCanopyEvaporation_dTCanair )%dat(1) ,& ! intent(in): [dp] derivative in canopy evaporation w.r.t. canopy air temperature - dCanopyEvaporation_dTCanopy => deriv_data%var(iLookDERIV%dCanopyEvaporation_dTCanopy )%dat(1) ,& ! intent(in): [dp] derivative in canopy evaporation w.r.t. canopy temperature - dCanopyEvaporation_dTGround => deriv_data%var(iLookDERIV%dCanopyEvaporation_dTGround )%dat(1) ,& ! intent(in): [dp] derivative in canopy evaporation w.r.t. ground temperature - dCanopyEvaporation_dCanWat => deriv_data%var(iLookDERIV%dCanopyEvaporation_dCanWat )%dat(1) ,& ! intent(in): [dp] derivative in canopy evaporation w.r.t. canopy total water content - dGroundEvaporation_dTCanair => deriv_data%var(iLookDERIV%dGroundEvaporation_dTCanair )%dat(1) ,& ! intent(in): [dp] derivative in ground evaporation w.r.t. canopy air temperature - dGroundEvaporation_dTCanopy => deriv_data%var(iLookDERIV%dGroundEvaporation_dTCanopy )%dat(1) ,& ! intent(in): [dp] derivative in ground evaporation w.r.t. canopy temperature - dGroundEvaporation_dTGround => deriv_data%var(iLookDERIV%dGroundEvaporation_dTGround )%dat(1) ,& ! intent(in): [dp] derivative in ground evaporation w.r.t. ground temperature - dGroundEvaporation_dCanWat => deriv_data%var(iLookDERIV%dGroundEvaporation_dCanWat )%dat(1) ,& ! intent(in): [dp] derivative in ground evaporation w.r.t. canopy total water content - ! derivatives in canopy water w.r.t canopy temperature - dCanLiq_dTcanopy => deriv_data%var(iLookDERIV%dCanLiq_dTcanopy )%dat(1) ,& ! intent(in): [dp] derivative in canopy liquid storage w.r.t. temperature - dTheta_dTkCanopy => deriv_data%var(iLookDERIV%dTheta_dTkCanopy )%dat(1) ,& ! intent(in): [dp] derivative in volumetric liquid water content w.r.t. temperature - d2Theta_dTkCanopy2 => deriv_data%var(iLookDERIV%d2Theta_dTkCanopy2 )%dat(1) ,& ! intent(in): [dp] second derivative of volumetric liquid water content w.r.t. temperature - dFracLiqVeg_dTkCanopy => deriv_data%var(iLookDERIV%dFracLiqVeg_dTkCanopy )%dat(1) ,& ! intent(in): [dp] derivative in fraction of (throughfall + drainage) w.r.t. temperature - ! derivatives in canopy liquid fluxes w.r.t. canopy water - scalarCanopyLiqDeriv => deriv_data%var(iLookDERIV%scalarCanopyLiqDeriv )%dat(1) ,& ! intent(in): [dp] derivative in (throughfall + drainage) w.r.t. canopy liquid water - ! derivatives in energy fluxes at the interface of snow+soil layers w.r.t. temperature in layers above and below - dNrgFlux_dTempAbove => deriv_data%var(iLookDERIV%dNrgFlux_dTempAbove )%dat ,& ! intent(in): [dp(:)] derivatives in the flux w.r.t. temperature in the layer above - dNrgFlux_dTempBelow => deriv_data%var(iLookDERIV%dNrgFlux_dTempBelow )%dat ,& ! intent(in): [dp(:)] derivatives in the flux w.r.t. temperature in the layer below - ! derivatives in energy fluxes at the interface of snow+soil layers w.r.t. water state in layers above and below - dNrgFlux_dWatAbove => deriv_data%var(iLookDERIV%dNrgFlux_dWatAbove )%dat ,& ! intent(in): [dp(:)] derivatives in the flux w.r.t. water state in the layer above - dNrgFlux_dWatBelow => deriv_data%var(iLookDERIV%dNrgFlux_dWatBelow )%dat ,& ! intent(in): [dp(:)] derivatives in the flux w.r.t. water state in the layer below - ! derivatives in soil transpiration w.r.t. canopy state variables - mLayerdTrans_dTCanair => deriv_data%var(iLookDERIV%mLayerdTrans_dTCanair )%dat ,& ! intent(in): derivatives in the soil layer transpiration flux w.r.t. canopy air temperature - mLayerdTrans_dTCanopy => deriv_data%var(iLookDERIV%mLayerdTrans_dTCanopy )%dat ,& ! intent(in): derivatives in the soil layer transpiration flux w.r.t. canopy temperature - mLayerdTrans_dTGround => deriv_data%var(iLookDERIV%mLayerdTrans_dTGround )%dat ,& ! intent(in): derivatives in the soil layer transpiration flux w.r.t. ground temperature - mLayerdTrans_dCanWat => deriv_data%var(iLookDERIV%mLayerdTrans_dCanWat )%dat ,& ! intent(in): derivatives in the soil layer transpiration flux w.r.t. canopy total water - ! derivatives in aquifer transpiration w.r.t. canopy state variables - dAquiferTrans_dTCanair => deriv_data%var(iLookDERIV%dAquiferTrans_dTCanair )%dat(1) ,& !intent(out): derivatives in the aquifer transpiration flux w.r.t. canopy air temperature - dAquiferTrans_dTCanopy => deriv_data%var(iLookDERIV%dAquiferTrans_dTCanopy )%dat(1) ,& ! intent(out): derivatives in the aquifer transpiration flux w.r.t. canopy temperature - dAquiferTrans_dTGround => deriv_data%var(iLookDERIV%dAquiferTrans_dTGround )%dat(1) ,& ! intent(out): derivatives in the aquifer transpiration flux w.r.t. ground temperature - dAquiferTrans_dCanWat => deriv_data%var(iLookDERIV%dAquiferTrans_dCanWat )%dat(1) ,& ! intent(out): derivatives in the aquifer transpiration flux w.r.t. canopy total water - ! derivative in liquid water fluxes at the interface of snow layers w.r.t. volumetric liquid water content in the layer above - iLayerLiqFluxSnowDeriv => deriv_data%var(iLookDERIV%iLayerLiqFluxSnowDeriv )%dat ,& ! intent(in): [dp(:)] derivative in vertical liquid water flux at layer interfaces - ! derivative in liquid water fluxes for the soil domain w.r.t hydrology state variables - dVolTot_dPsi0 => deriv_data%var(iLookDERIV%dVolTot_dPsi0 )%dat ,& ! intent(in): [dp(:)] derivative in total water content w.r.t. total water matric potential - d2VolTot_d2Psi0 => deriv_data%var(iLookDERIV%d2VolTot_d2Psi0 )%dat ,& ! intent(in): [dp(:)] second derivative in total water content w.r.t. total water matric potential - dCompress_dPsi => deriv_data%var(iLookDERIV%dCompress_dPsi )%dat ,& ! intent(in): [dp(:)] derivative in compressibility w.r.t matric head - dq_dHydStateAbove => deriv_data%var(iLookDERIV%dq_dHydStateAbove )%dat ,& ! intent(in): [dp(:)] change in flux at layer interfaces w.r.t. states in the layer above - dq_dHydStateBelow => deriv_data%var(iLookDERIV%dq_dHydStateBelow )%dat ,& ! intent(in): [dp(:)] change in flux at layer interfaces w.r.t. states in the layer below - dq_dHydStateLayerSurfVec => deriv_data%var(iLookDERIV%dq_dHydStateLayerSurfVec )%dat ,& ! intent(in): [dp(:)] change in the flux in soil surface interface w.r.t. state variables in layers - ! derivative in baseflow flux w.r.t. aquifer storage - dBaseflow_dAquifer => deriv_data%var(iLookDERIV%dBaseflow_dAquifer )%dat(1) ,& ! intent(out): [dp(:)] derivative in baseflow flux w.r.t. aquifer storage (s-1) - ! derivative in liquid water fluxes for the soil domain w.r.t energy state variables - dq_dNrgStateAbove => deriv_data%var(iLookDERIV%dq_dNrgStateAbove )%dat ,& ! intent(in): [dp(:)] change in flux at layer interfaces w.r.t. states in the layer above - dq_dNrgStateBelow => deriv_data%var(iLookDERIV%dq_dNrgStateBelow )%dat ,& ! intent(in): [dp(:)] change in flux at layer interfaces w.r.t. states in the layer below - dq_dNrgStateLayerSurfVec => deriv_data%var(iLookDERIV%dq_dNrgStateLayerSurfVec )%dat ,& ! intent(in): [dp(:)] change in the flux in soil surface interface w.r.t. state variables in layers - ! derivative in liquid water fluxes for the soil and snow domain w.r.t temperature - dFracLiqSnow_dTk => deriv_data%var(iLookDERIV%dFracLiqSnow_dTk )%dat ,& ! intent(in): [dp(:)] derivative in fraction of liquid snow w.r.t. temperature - mLayerdTheta_dTk => deriv_data%var(iLookDERIV%mLayerdTheta_dTk )%dat ,& ! intent(in): [dp(:)] derivative in volumetric liquid water content w.r.t. temperature - mLayerd2Theta_dTk2 => deriv_data%var(iLookDERIV%mLayerd2Theta_dTk2 )%dat ,& ! intent(in): [dp(:)] second derivative of volumetric liquid water content w.r.t. temperature - ! derivate in bulk heat capacity w.r.t. relevant state variables - dVolHtCapBulk_dPsi0 => deriv_data%var(iLookDERIV%dVolHtCapBulk_dPsi0 )%dat ,& ! intent(in): [dp(:)] derivative in bulk heat capacity w.r.t. matric potential - dVolHtCapBulk_dTheta => deriv_data%var(iLookDERIV%dVolHtCapBulk_dTheta )%dat ,& ! intent(in): [dp(:)] derivative in bulk heat capacity w.r.t. volumetric water content - dVolHtCapBulk_dCanWat => deriv_data%var(iLookDERIV%dVolHtCapBulk_dCanWat )%dat(1) ,& ! intent(in): [dp] derivative in bulk heat capacity w.r.t. volumetric water content - dVolHtCapBulk_dTk => deriv_data%var(iLookDERIV%dVolHtCapBulk_dTk )%dat ,& ! intent(in): [dp(:)] derivative in bulk heat capacity w.r.t. temperature - dVolHtCapBulk_dTkCanopy => deriv_data%var(iLookDERIV%dVolHtCapBulk_dTkCanopy )%dat(1) ,& ! intent(in): [dp] derivative in bulk heat capacity w.r.t. temperature - ! diagnostic variables - scalarFracLiqVeg => diag_data%var(iLookDIAG%scalarFracLiqVeg )%dat(1) ,& ! intent(in): [dp] fraction of liquid water on vegetation (-) - scalarBulkVolHeatCapVeg => diag_data%var(iLookDIAG%scalarBulkVolHeatCapVeg )%dat(1) ,& ! intent(in): [dp] bulk volumetric heat capacity of vegetation (J m-3 K-1) - mLayerFracLiqSnow => diag_data%var(iLookDIAG%mLayerFracLiqSnow )%dat ,& ! intent(in): [dp(:)] fraction of liquid water in each snow layer (-) - mLayerVolHtCapBulk => diag_data%var(iLookDIAG%mLayerVolHtCapBulk )%dat ,& ! intent(in): [dp(:)] bulk volumetric heat capacity in each snow and soil layer (J m-3 K-1) - scalarSoilControl => diag_data%var(iLookDIAG%scalarSoilControl )%dat(1) ,& ! intent(in): [dp] soil control on infiltration, zero or one - ! canopy and layer depth - canopyDepth => diag_data%var(iLookDIAG%scalarCanopyDepth )%dat(1) ,& ! intent(in): [dp ] canopy depth (m) - mLayerDepth => prog_data%var(iLookPROG%mLayerDepth )%dat ,& ! intent(in): [dp(:)] depth of each layer in the snow-soil sub-domain (m) - layerType => indx_data%var(iLookINDEX%layerType )%dat & ! intent(in): [i4b(:)] named variables defining the type of layer in snow+soil domain - ) ! making association with data in structures - ! -------------------------------------------------------------- - ! initialize error control - err=0; message='computJacDAE/' - - ! ********************************************************************************************************************************************************* - ! ********************************************************************************************************************************************************* - ! * PART 0: PRELIMINARIES (INITIALIZE JACOBIAN AND COMPUTE TIME-VARIABLE DIAGONAL TERMS) - ! ********************************************************************************************************************************************************* - ! ********************************************************************************************************************************************************* - - ! get the number of state variables - nState = size(dMat) - - ! initialize the Jacobian - ! NOTE: this needs to be done every time, since Jacobian matrix is modified in the solver - aJac(:,:) = 0._rkind ! analytical Jacobian matrix - - ! compute terms in the Jacobian for vegetation (excluding fluxes) - ! NOTE: energy for vegetation is computed *within* the iteration loop as it includes phase change - if(ixVegNrg/=integerMissing)then - dMat(ixVegNrg) = ( scalarBulkVolHeatCapVeg + LH_fus*iden_water*dTheta_dTkCanopy ) * cj & - + dVolHtCapBulk_dTkCanopy * scalarCanopyTempPrime & - + LH_fus*iden_water * scalarCanopyTempPrime * d2Theta_dTkCanopy2 & - + LH_fus * dFracLiqVeg_dTkCanopy * scalarCanopyWatPrime / canopyDepth - end if - - ! compute additional terms for the Jacobian for the snow-soil domain (excluding fluxes) - ! NOTE: energy for snow+soil is computed *within* the iteration loop as it includes phase change - do iLayer=1,nLayers - if(ixSnowSoilNrg(iLayer)/=integerMissing) then - dMat(ixSnowSoilNrg(iLayer)) = ( mLayerVolHtCapBulk(iLayer) + LH_fus*iden_water*mLayerdTheta_dTk(iLayer) ) * cj & - + dVolHtCapBulk_dTk(iLayer) * mLayerTempPrime(iLayer) & - + LH_fus*iden_water * mLayerTempPrime(iLayer) * mLayerd2Theta_dTk2(iLayer) & - + LH_fus*iden_water * dFracLiqSnow_dTk(iLayer) * mLayerVolFracWatPrime(iLayer) - end if - end do - - ! compute additional terms for the Jacobian for the soil domain (excluding fluxes) - do iLayer=1,nSoil - if(ixSoilOnlyHyd(iLayer)/=integerMissing)then - dMat(ixSoilOnlyHyd(iLayer)) = ( dVolTot_dPsi0(iLayer) + dCompress_dPsi(iLayer) ) * cj + d2VolTot_d2Psi0(iLayer) * mLayerMatricHeadPrime(iLayer) - - if(ixRichards==mixdform)then - dMat(ixSoilOnlyHyd(iLayer)) = dMat(ixSoilOnlyHyd(iLayer)) + specificStorage * dVolTot_dPsi0(iLayer) * mLayerMatricHeadPrime(iLayer) / theta_sat(iLayer) - end if - - end if - end do - - ! define the form of the matrix - select case(ixMatrix) - - ! ********************************************************************************************************************************************************* - ! ********************************************************************************************************************************************************* - ! * PART 1: BAND MATRIX - ! ********************************************************************************************************************************************************* - ! ********************************************************************************************************************************************************* - case(ixBandMatrix) ! ixBandMatrix ixFullMatrix - print *, 'banded jacobian matrix needs to be implemented' - stop 1 - - ! ********************************************************************************************************************************************************* - ! ********************************************************************************************************************************************************* - ! * PART 2: FULL MATRIX - ! ********************************************************************************************************************************************************* - ! ********************************************************************************************************************************************************* - case(ixFullMatrix) ! ixFullMatrix ixBandMatrix - - ! check - if(size(aJac,1)/=size(dMat) .or. size(aJac,2)/=size(dMat))then - message=trim(message)//'unexpected shape of the Jacobian matrix: expect aJac(nState,nState)' - err=20; return - end if - - ! ----- - ! * energy and liquid fluxes over vegetation... - ! --------------------------------------------- - if(computeVegFlux)then ! (derivatives only defined when vegetation protrudes over the surface) - - ! * liquid water fluxes for vegetation canopy (-), dt*scalarFracLiqVeg*scalarCanopyLiqDeriv is the derivative in throughfall and canopy drainage with canopy water - if(ixVegHyd/=integerMissing) aJac(ixVegHyd,ixVegHyd) = -scalarFracLiqVeg*(dCanopyEvaporation_dCanWat - scalarCanopyLiqDeriv)*dt + 1._rkind * cj - - ! * cross-derivative terms for canopy water - if(ixVegHyd/=integerMissing)then - - ! cross-derivative terms w.r.t. system temperatures (kg m-2 K-1) - if(ixCasNrg/=integerMissing) aJac(ixVegHyd,ixCasNrg) = -dCanopyEvaporation_dTCanair*dt - ! dt*scalarCanopyLiqDeriv*dCanLiq_dTcanopy is the derivative in throughfall and canopy drainage with canopy temperature - if(ixVegNrg/=integerMissing) aJac(ixVegHyd,ixVegNrg) = -dCanopyEvaporation_dTCanopy*dt + dt*scalarCanopyLiqDeriv*dCanLiq_dTcanopy - if(ixTopNrg/=integerMissing) aJac(ixVegHyd,ixTopNrg) = -dCanopyEvaporation_dTGround*dt - - ! cross-derivative terms w.r.t. canopy water (kg-1 m2) - if(ixTopHyd/=integerMissing) aJac(ixTopHyd,ixVegHyd) = (dt/mLayerDepth(1))*(-scalarSoilControl*scalarFracLiqVeg*scalarCanopyLiqDeriv)/iden_water - - ! cross-derivative terms w.r.t. canopy liquid water (J m-1 kg-1) - ! NOTE: dIce/dLiq = (1 - scalarFracLiqVeg); dIce*LH_fus/canopyDepth = J m-3; dLiq = kg m-2 - if(ixVegNrg/=integerMissing) aJac(ixVegNrg,ixVegHyd) = (-1._rkind + scalarFracLiqVeg)*LH_fus/canopyDepth * cj & - + dVolHtCapBulk_dCanWat * scalarCanopyTempPrime - (dt/canopyDepth) * dCanopyNetFlux_dCanWat & - + LH_fus * scalarCanopyTempPrime * dFracLiqVeg_dTkCanopy / canopyDepth ! dF/dLiq - if(ixTopNrg/=integerMissing) aJac(ixTopNrg,ixVegHyd) = (dt/mLayerDepth(1))*(-dGroundNetFlux_dCanWat) - endif - - ! cross-derivative terms w.r.t. canopy temperature (K-1) - if(ixVegNrg/=integerMissing)then - if(ixTopHyd/=integerMissing) aJac(ixTopHyd,ixVegNrg) = (dt/mLayerDepth(1))*(-scalarSoilControl*scalarCanopyLiqDeriv*dCanLiq_dTcanopy)/iden_water - endif - - ! energy fluxes with the canopy air space (J m-3 K-1) - if(ixCasNrg/=integerMissing)then - aJac(ixCasNrg,ixCasNrg) = (dt/canopyDepth)*(-dCanairNetFlux_dCanairTemp) + dMat(ixCasNrg) * cj - if(ixVegNrg/=integerMissing) aJac(ixCasNrg,ixVegNrg) = (dt/canopyDepth)*(-dCanairNetFlux_dCanopyTemp) - if(ixTopNrg/=integerMissing) aJac(ixCasNrg,ixTopNrg) = (dt/canopyDepth)*(-dCanairNetFlux_dGroundTemp) - endif - - ! energy fluxes with the vegetation canopy (J m-3 K-1) - if(ixVegNrg/=integerMissing)then - if(ixCasNrg/=integerMissing) aJac(ixVegNrg,ixCasNrg) = (dt/canopyDepth)*(-dCanopyNetFlux_dCanairTemp) - aJac(ixVegNrg,ixVegNrg) = (dt/canopyDepth)*(-dCanopyNetFlux_dCanopyTemp) + dMat(ixVegNrg) - if(ixTopNrg/=integerMissing) aJac(ixVegNrg,ixTopNrg) = (dt/canopyDepth)*(-dCanopyNetFlux_dGroundTemp) - endif - - ! energy fluxes with the surface (J m-3 K-1) - if(ixTopNrg/=integerMissing)then - if(ixCasNrg/=integerMissing) aJac(ixTopNrg,ixCasNrg) = (dt/mLayerDepth(1))*(-dGroundNetFlux_dCanairTemp) - if(ixVegNrg/=integerMissing) aJac(ixTopNrg,ixVegNrg) = (dt/mLayerDepth(1))*(-dGroundNetFlux_dCanopyTemp) - endif - - endif ! if there is a need to compute energy fluxes within vegetation - - ! ----- - ! * energy fluxes for the snow+soil domain... - ! ------------------------------------------- - if(nSnowSoilNrg>0)then - do iLayer=1,nLayers ! loop through all layers in the snow+soil domain - - ! check if the state is in the subset - if(ixSnowSoilNrg(iLayer)==integerMissing) cycle - - ! - define index within the state subset and the full state vector - jState = ixSnowSoilNrg(iLayer) ! index within the state subset - - ! - diagonal elements - aJac(jState,jState) = (dt/mLayerDepth(iLayer))*(-dNrgFlux_dTempBelow(iLayer-1) + dNrgFlux_dTempAbove(iLayer)) + dMat(jState) - - ! - lower-diagonal elements - if(iLayer > 1)then - if(ixSnowSoilNrg(iLayer-1)/=integerMissing) aJac(ixSnowSoilNrg(iLayer-1),jState) = (dt/mLayerDepth(iLayer-1))*( dNrgFlux_dTempBelow(iLayer-1) ) - endif - - ! - upper diagonal elements - if(iLayer < nLayers)then - if(ixSnowSoilNrg(iLayer+1)/=integerMissing) aJac(ixSnowSoilNrg(iLayer+1),jState) = (dt/mLayerDepth(iLayer+1))*(-dNrgFlux_dTempAbove(iLayer ) ) - endif - - end do ! (looping through energy states in the snow+soil domain) - endif ! (if the subset includes energy state variables in the snow+soil domain) - - ! ----- - ! * liquid water fluxes for the snow domain... - ! -------------------------------------------- - if(nSnowOnlyHyd>0)then - do iLayer=1,nSnow ! loop through layers in the snow domain - - ! - check that the snow layer is desired - if(ixSnowOnlyHyd(iLayer)==integerMissing) cycle - - ! - define state indices for the current layer - watState = ixSnowOnlyHyd(iLayer) ! hydrology state index within the state subset - - ! compute factor to convert liquid water derivative to total water derivative - select case( ixHydType(iLayer) ) - case(iname_watLayer); convLiq2tot = mLayerFracLiqSnow(iLayer) - case default; convLiq2tot = 1._rkind - end select - - ! - diagonal elements - aJac(watState,watState) = (dt/mLayerDepth(iLayer))*iLayerLiqFluxSnowDeriv(iLayer)*convLiq2tot + dMat(watState) * cj - - ! - lower-diagonal elements - if(iLayer > 1)then - if(ixSnowOnlyHyd(iLayer-1)/=integerMissing) aJac(ixSnowOnlyHyd(iLayer-1),watState) = 0._rkind ! sub-diagonal: no dependence on other layers - endif - - ! - upper diagonal elements - if(iLayer < nSnow)then - if(ixSnowOnlyHyd(iLayer+1)/=integerMissing) aJac(ixSnowOnlyHyd(iLayer+1),watState) = -(dt/mLayerDepth(iLayer+1))*iLayerLiqFluxSnowDeriv(iLayer)*convLiq2tot ! dVol(below)/dLiq(above) -- (-) - endif - - end do ! (looping through liquid water states in the snow domain) - endif ! (if the subset includes hydrology state variables in the snow domain) - - ! ----- - ! * cross derivatives in the snow domain... - ! ---------------------------------------- - if(nSnowOnlyHyd>0 .and. nSnowOnlyNrg>0)then - do iLayer=1,nSnow ! loop through layers in the snow domain - - ! - check that the snow layer is desired - if(ixSnowOnlyNrg(iLayer)==integerMissing) cycle - - ! (define the energy state) - nrgState = ixSnowOnlyNrg(iLayer) ! index within the full state vector - - ! - define state indices for the current layer - watState = ixSnowOnlyHyd(iLayer) ! hydrology state index within the state subset - - if(watstate/=integerMissing)then ! (energy state for the current layer is within the state subset) - - ! - include derivatives of energy fluxes w.r.t water fluxes for current layer - aJac(nrgState,watState) = (-1._rkind + mLayerFracLiqSnow(iLayer))*LH_fus*iden_water * cj & - + dVolHtCapBulk_dTheta(iLayer) * mLayerTempPrime(iLayer) & - + (dt/mLayerDepth(iLayer))*(-dNrgFlux_dWatBelow(iLayer-1) + dNrgFlux_dWatAbove(iLayer)) & - + LH_fus*iden_water * mLayerTempPrime(iLayer) * dFracLiqSnow_dTk(iLayer) ! (dF/dLiq) - - ! - include derivatives of water fluxes w.r.t energy fluxes for current layer - aJac(watState,nrgState) = (dt/mLayerDepth(iLayer))*iLayerLiqFluxSnowDeriv(iLayer)*mLayerdTheta_dTk(iLayer) ! (dVol/dT) - - ! (cross-derivative terms for the layer below) - if(iLayer<nSnow)then - aJac(ixSnowOnlyHyd(iLayer+1),nrgState) = -(dt/mLayerDepth(iLayer+1))*iLayerLiqFluxSnowDeriv(iLayer)*mLayerdTheta_dTk(iLayer) ! dVol(below)/dT(above) -- K-1 - endif ! (if there is a water state in the layer below the current layer in the given state subset) - - ! - include derivatives of heat capacity w.r.t water fluxes for surrounding layers starting with layer above - if(iLayer>1)then - if(ixSnowOnlyNrg(iLayer-1)/=integerMissing) aJac(ixSnowOnlyNrg(iLayer-1),watState) = (dt/mLayerDepth(iLayer-1))*( dNrgFlux_dWatBelow(iLayer-1) ) - endif - - ! (cross-derivative terms for the layer below) - if(iLayer<nSnow)then - if(ixSnowOnlyNrg(iLayer+1)/=integerMissing) aJac(ixSnowOnlyNrg(iLayer+1),watState) = (dt/mLayerDepth(iLayer+1))*(-dNrgFlux_dWatAbove(iLayer ) ) - elseif(iLayer==nSnow .and. nSoilOnlyNrg>0)then !bottom snow layer and there is soil below - if(ixSoilOnlyNrg(1)/=integerMissing) aJac(ixSoilOnlyNrg(1),watState) = (dt/mLayerDepth(nSnow+1))*(-dNrgFlux_dWatAbove(nSnow) ) - endif - - endif ! (if the energy state for the current layer is within the state subset) - - end do ! (looping through snow layers) - - endif ! (if there are state variables for both water and energy in the snow domain) - - ! ----- - ! * liquid water fluxes for the soil domain... - ! -------------------------------------------- - if(nSoilOnlyHyd>0)then - do iLayer=1,nSoil - - ! - check that the soil layer is desired - if(ixSoilOnlyHyd(iLayer)==integerMissing) cycle - - ! - define state indices - watState = ixSoilOnlyHyd(iLayer) ! hydrology state index within the state subset - - ! - define indices of the soil layers - jLayer = iLayer+nSnow ! index of layer in the snow+soil vector - - ! - compute the diagonal elements - ! all terms *excluding* baseflow - aJac(watState,watState) = (dt/mLayerDepth(jLayer))*(-dq_dHydStateBelow(iLayer-1) + dq_dHydStateAbove(iLayer)) + dMat(watState) - - ! - compute the lower-diagonal elements - if(iLayer > 1)then - if(ixSoilOnlyHyd(iLayer-1)/=integerMissing) aJac(ixSoilOnlyHyd(iLayer-1),watState) = (dt/mLayerDepth(jLayer-1))*( dq_dHydStateBelow(iLayer-1)) - endif - - ! - compute the upper-diagonal elements - if(iLayer<nSoil)then - if(ixSoilOnlyHyd(iLayer+1)/=integerMissing) aJac(ixSoilOnlyHyd(iLayer+1),watState) = (dt/mLayerDepth(jLayer+1))*(-dq_dHydStateAbove(iLayer)) - endif - - ! - include terms for baseflow - if(computeBaseflow .and. nSoilOnlyHyd==nSoil)then - do pLayer=1,nSoil - qState = ixSoilOnlyHyd(pLayer) ! hydrology state index within the state subset - aJac(watState,qState) = (dt/mLayerDepth(jLayer))*dBaseflow_dMatric(iLayer,pLayer) + aJac(watState,qState) - end do - endif - - ! - include terms for surface infiltration below surface - if(ixSoilOnlyHyd(1)/=integerMissing) aJac(ixSoilOnlyHyd(1),watState) = -(dt/mLayerDepth(1+nSnow))*dq_dHydStateLayerSurfVec(iLayer) + aJac(ixSoilOnlyHyd(1),watState) - - end do ! (looping through hydrology states in the soil domain) - - ! - include terms for surface infiltration above surface - if(nSnowOnlyHyd>0 .and. ixSnowOnlyHyd(nSnow)/=integerMissing)then - if(ixSoilOnlyHyd(1)/=integerMissing) aJac(ixSoilOnlyHyd(1),ixSnowOnlyHyd(nSnow)) = -(dt/mLayerDepth(1+nSnow))*dq_dHydStateLayerSurfVec(0) - elseif(computeVegFlux .and. ixVegHyd/=integerMissing)then !ixTopHyd = ixSoilOnlyHyd(1) - if(ixTopHyd/=integerMissing) aJac(ixTopHyd,ixVegHyd) = -(dt/mLayerDepth(1+nSnow))*dq_dHydStateLayerSurfVec(0) + aJac(ixTopHyd,ixVegNrg) - endif - - endif ! (if the subset includes hydrology state variables in the soil domain) - - ! ----- - ! * liquid water fluxes for the aquifer... - ! ---------------------------------------- - if(ixAqWat/=integerMissing) then - aJac(ixAqWat,ixAqWat) = -dBaseflow_dAquifer*dt + dMat(ixAqWat) * cj - aJac(ixAqWat,ixSoilOnlyHyd(nSoil)) = -dq_dHydStateAbove(nSoil)*dt ! dAquiferRecharge_dWat = d_iLayerLiqFluxSoil(nSoil)_dWat - aJac(ixAqWat,ixSoilOnlyNrg(nSoil)) = -dq_dNrgStateAbove(nSoil)*dt ! dAquiferRecharge_dTk = d_iLayerLiqFluxSoil(nSoil)_dTk - ! - include derivatives of energy and water w.r.t soil transpiration (dependent on canopy transpiration) - if(computeVegFlux)then - aJac(ixAqWat,ixVegHyd) = -dAquiferTrans_dCanWat*dt ! dVol/dLiq (kg m-2)-1 - aJac(ixAqWat,ixCasNrg) = -dAquiferTrans_dTCanair*dt ! dVol/dT (K-1) - aJac(ixAqWat,ixVegNrg) = -dAquiferTrans_dTCanopy*dt ! dVol/dT (K-1) - aJac(ixAqWat,ixTopNrg) = -dAquiferTrans_dTGround*dt ! dVol/dT (K-1) - endif - endif - - ! ----- - ! * cross derivatives in the soil domain... - ! ---------------------------------------- - if(nSoilOnlyHyd>0 .and. nSoilOnlyNrg>0)then - do iLayer=1,nSoilOnlyNrg - - ! - check that the soil layer is desired - if(ixSoilOnlyNrg(iLayer)==integerMissing) cycle - - ! - define indices of the soil layers - jLayer = iLayer+nSnow ! index of layer in the snow+soil vector - - ! - define the energy state variable - nrgState = ixNrgLayer(jLayer) ! index within the full state vector - - ! - define index of hydrology state variable within the state subset - watState = ixSoilOnlyHyd(iLayer) - - ! only compute derivatives if the energy state for the current layer is within the state subset - if(watstate/=integerMissing)then - - ! - include derivatives in liquid water fluxes w.r.t. temperature for current layer - aJac(watState,nrgState) = (dt/mLayerDepth(jLayer))*(-dq_dNrgStateBelow(iLayer-1) + dq_dNrgStateAbove(iLayer)) ! dVol/dT (K-1) -- flux depends on ice impedance - - ! - compute lower diagonal elements - if(iLayer>1)then - if(ixSoilOnlyHyd(iLayer-1)/=integerMissing) aJac(ixSoilOnlyHyd(iLayer-1),nrgState) = (dt/mLayerDepth(jLayer-1))*( dq_dNrgStateBelow(iLayer-1)) ! K-1 - endif - - ! compute upper-diagonal elements - if(iLayer<nSoil)then - if(ixSoilOnlyHyd(iLayer+1)/=integerMissing) aJac(ixSoilOnlyHyd(iLayer+1),nrgState) = (dt/mLayerDepth(jLayer+1))*(-dq_dNrgStateAbove(iLayer)) ! K-1 - endif - - ! - include derivatives of energy w.r.t. ground evaporation - if(nSnow==0 .and. iLayer==1)then ! upper-most soil layer - if(computeVegFlux)then - aJac(watState,ixVegHyd) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dCanWat/iden_water) ! dVol/dLiq (kg m-2)-1 - aJac(watState,ixCasNrg) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dTCanair/iden_water) ! dVol/dT (K-1) - aJac(watState,ixVegNrg) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dTCanopy/iden_water) ! dVol/dT (K-1) - endif - aJac(watState,ixTopNrg) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dTGround/iden_water) + aJac(watState,ixTopNrg) ! dVol/dT (K-1) - endif - - ! - include derivatives of energy and water w.r.t soil transpiration (dependent on canopy transpiration) - if(computeVegFlux)then - aJac(watState,ixVegHyd) = (dt/mLayerDepth(jLayer))*(-mLayerdTrans_dCanWat(iLayer)) + aJac(watState,ixVegHyd) ! dVol/dLiq (kg m-2)-1 - aJac(watState,ixCasNrg) = (dt/mLayerDepth(jLayer))*(-mLayerdTrans_dTCanair(iLayer)) + aJac(watState,ixCasNrg) ! dVol/dT (K-1) - aJac(watState,ixVegNrg) = (dt/mLayerDepth(jLayer))*(-mLayerdTrans_dTCanopy(iLayer)) + aJac(watState,ixVegNrg) ! dVol/dT (K-1) - aJac(watState,ixTopNrg) = (dt/mLayerDepth(jLayer))*(-mLayerdTrans_dTGround(iLayer)) + aJac(watState,ixTopNrg) ! dVol/dT (K-1) - endif - - ! - include derivatives in energy fluxes w.r.t. with respect to water for current layer - aJac(nrgState,watState) = dVolHtCapBulk_dPsi0(iLayer) * mLayerTempPrime(jLayer) & - + (dt/mLayerDepth(jLayer))*(-dNrgFlux_dWatBelow(jLayer-1) + dNrgFlux_dWatAbove(jLayer)) - if(mLayerdTheta_dTk(jLayer) > verySmall)then ! ice is present - aJac(nrgState,watState) = -LH_fus*iden_water * dVolTot_dPsi0(iLayer) * cj & - - LH_fus*iden_water * mLayerMatricHeadPrime(iLayer) * d2VolTot_d2Psi0(iLayer) + aJac(nrgState,watState) ! dNrg/dMat (J m-3 m-1) -- dMat changes volumetric water, and hence ice content - endif - - ! - include derivatives of heat capacity w.r.t water fluxes for surrounding layers starting with layer above - if(iLayer>1)then - if(ixSoilOnlyNrg(iLayer-1)/=integerMissing) aJac(ixSoilOnlyNrg(iLayer-1),watState) = (dt/mLayerDepth(jLayer-1))*( dNrgFlux_dWatBelow(jLayer-1) ) - elseif(iLayer==1 .and. nSnowOnlyNrg>0)then !top soil layer and there is snow above - if(ixSnowOnlyNrg(nSnow)/=integerMissing) aJac(ixSnowOnlyNrg(nSnow),watState) = (dt/mLayerDepth(nSnow))*( dNrgFlux_dWatBelow(nSnow) ) - endif - - ! (cross-derivative terms for the layer below) - if(iLayer<nSoil)then - if(ixSoilOnlyHyd(iLayer+1)/=integerMissing) aJac(ixSoilOnlyNrg(iLayer+1),watState) = (dt/mLayerDepth(jLayer+1))*(-dNrgFlux_dWatAbove(jLayer ) ) - endif - - endif ! (if the water state for the current layer is within the state subset) - - ! - include terms for surface infiltration below surface - if(ixSoilOnlyHyd(1)/=integerMissing) aJac(ixSoilOnlyHyd(1),nrgState) = -(dt/mLayerDepth(1+nSnow))*dq_dNrgStateLayerSurfVec(iLayer) + aJac(ixSoilOnlyHyd(1),nrgState) - - end do ! (looping through soil layers) - - ! - include terms for surface infiltration above surface - if(nSnowOnlyHyd>0 .and. ixSnowOnlyNrg(nSnow)/=integerMissing)then - if(ixSoilOnlyHyd(1)/=integerMissing) aJac(ixSoilOnlyHyd(1),ixSnowOnlyNrg(nSnow)) = -(dt/mLayerDepth(1+nSnow))*dq_dNrgStateLayerSurfVec(0) - elseif(computeVegFlux .and. ixVegNrg/=integerMissing) then !ixTopHyd = ixSoilOnlyHyd(1) - if(ixTopHyd/=integerMissing) aJac(ixTopHyd,ixVegNrg) = -(dt/mLayerDepth(1+nSnow))*dq_dNrgStateLayerSurfVec(0) + aJac(ixTopHyd,ixVegNrg) - endif - - endif ! (if there are state variables for both water and energy in the soil domain) - - ! print the Jacobian - if(globalPrintFlag)then - print*, '** analytical Jacobian (full):' - write(*,'(a4,1x,100(i12,1x))') 'xCol', (iLayer, iLayer=min(iJac1,nState),min(iJac2,nState)) - do iLayer=min(iJac1,nState),min(iJac2,nState) - write(*,'(i4,1x,100(e12.5,1x))') iLayer, aJac(min(iJac1,nState):min(iJac2,nState),iLayer) - end do - end if - - !print*, '** analytical Jacobian (full):' - !write(*,'(a4,1x,100(i12,1x))') 'xCol', (iLayer, iLayer=1,size(aJac,2)) - !do iLayer=1,size(aJac,1) - ! write(*,'(i4,1x,100(e12.5,1x))') iLayer, aJac(iLayer,1:size(aJac,2)) - !end do - !print *, '--------------------------------------------------------------' - - ! *** - ! check - case default; err=20; message=trim(message)//'unable to identify option for the type of matrix'; return - - end select ! type of matrix - - if(any(isNan(aJac)))then - print *, '******************************* WE FOUND NAN IN JACOBIAN ************************************' - stop 1 - message=trim(message)//'we found NaN' - err=20; return - endif - - ! end association to variables in the data structures - end associate - - end subroutine computJacDAE - - - ! ********************************************************************************************************** - ! private function: get the off-diagonal index in the band-diagonal matrix - ! ********************************************************************************************************** - function ixOffDiag(jState,iState) - implicit none - integer(i4b),intent(in) :: jState ! off-diagonal state - integer(i4b),intent(in) :: iState ! diagonal state - integer(i4b) :: ixOffDiag ! off-diagonal index in gthe band-diagonal matrix - ixOffDiag = ixBandOffset + jState - iState - end function ixOffDiag - -end module computJacDAE_module diff --git a/build/source/engine/sundials/computJacobSundials.f90 b/build/source/engine/sundials/computJacobSundials.f90 new file mode 100644 index 0000000..e0f6df8 --- /dev/null +++ b/build/source/engine/sundials/computJacobSundials.f90 @@ -0,0 +1,1120 @@ +! SUMMA - Structure for Unifying Multiple Modeling Alternatives +! Copyright (C) 2014-2015 NCAR/RAL +! +! This file is part of SUMMA +! +! For more information see: http://www.ral.ucar.edu/projects/summa +! +! This program is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with this program. If not, see <http://www.gnu.org/licenses/>. + +module computJacobSundials_module + + ! data types + USE nrtype + USE type4IDA + + ! derived types to define the data structures + USE data_types,only:& + var_i, & ! data vector (i4b) + var_d, & ! data vector (rkind) + var_ilength, & ! data vector with variable length dimension (i4b) + var_dlength, & ! data vector with variable length dimension (rkind) + model_options ! defines the model decisions + + ! named variables for structure elements + USE var_lookup,only:iLookDECISIONS ! named variables for elements of the decision structure + USE var_lookup,only:iLookPARAM ! named variables for structure elements + USE var_lookup,only:iLookPROG ! named variables for structure elements + USE var_lookup,only:iLookINDEX ! named variables for structure elements + USE var_lookup,only:iLookDIAG ! named variables for structure elements + USE var_lookup,only:iLookFLUX ! named variables for structure elements + USE var_lookup,only:iLookDERIV ! named variables for structure elements + + ! access the global print flag + USE globalData,only:globalPrintFlag + + ! access missing values + USE globalData,only:integerMissing ! missing integer + USE globalData,only:realMissing ! missing real number + USE globalData,only:quadMissing ! missing quadruple precision number + + ! domain types + USE globalData,only:iname_veg ! named variables for vegetation + USE globalData,only:iname_snow ! named variables for snow + USE globalData,only:iname_soil ! named variables for soil + + ! named variables to describe the state variable type + USE globalData,only:iname_nrgCanair ! named variable defining the energy of the canopy air space + USE globalData,only:iname_nrgCanopy ! named variable defining the energy of the vegetation canopy + USE globalData,only:iname_watCanopy ! named variable defining the mass of water on the vegetation canopy + USE globalData,only:iname_nrgLayer ! named variable defining the energy state variable for snow+soil layers + USE globalData,only:iname_watLayer ! named variable defining the total water state variable for snow+soil layers + USE globalData,only:iname_liqLayer ! named variable defining the liquid water state variable for snow+soil layers + USE globalData,only:iname_matLayer ! named variable defining the matric head state variable for soil layers + USE globalData,only:iname_lmpLayer ! named variable defining the liquid matric potential state variable for soil layers + USE globalData,only:model_decisions ! model decision structure + + ! access named variables to describe the form and structure of the matrices used in the numerical solver + USE globalData,only: ku ! number of super-diagonal bands + USE globalData,only: kl ! number of sub-diagonal bands + USE globalData,only: ixDiag ! index for the diagonal band + USE globalData,only: nBands ! length of the leading dimension of the band diagonal matrix + USE globalData,only: ixFullMatrix ! named variable for the full Jacobian matrix + USE globalData,only: ixBandMatrix ! named variable for the band diagonal matrix + USE globalData,only: iJac1 ! first layer of the Jacobian to print + USE globalData,only: iJac2 ! last layer of the Jacobian to print + + ! constants + USE multiconst,only:& + Tfreeze, & ! temperature at freezing (K) + LH_fus, & ! latent heat of fusion (J kg-1) + LH_vap, & ! latent heat of vaporization (J kg-1) + LH_sub, & ! latent heat of sublimation (J kg-1) + Cp_air, & ! specific heat of air (J kg-1 K-1) + iden_air, & ! intrinsic density of air (kg m-3) + iden_ice, & ! intrinsic density of ice (kg m-3) + iden_water ! intrinsic density of liquid water (kg m-3) + + ! look-up values for the choice of groundwater representation (local-column, or single-basin) + USE mDecisions_module,only: & + localColumn, & ! separate groundwater representation in each local soil column + singleBasin ! single groundwater store over the entire basin + + ! look-up values for the choice of groundwater parameterization + USE mDecisions_module,only: & + qbaseTopmodel, & ! TOPMODEL-ish baseflow parameterization + bigBucket, & ! a big bucket (lumped aquifer model) + noExplicit ! no explicit groundwater parameterization + + ! look-up values for the form of Richards' equation + USE mDecisions_module,only: & + moisture, & ! moisture-based form of Richards' equation + mixdform ! mixed form of Richards' equation + + use, intrinsic :: iso_c_binding + + implicit none + ! define constants + real(rkind),parameter :: verySmall=tiny(1.0_rkind) ! a very small number + integer(i4b),parameter :: ixBandOffset=kl+ku+1 ! offset in the band Jacobian matrix + + private + public::computJacobSundials + public::computJacobSetup + public::computJacob4IDA + contains + + ! ********************************************************************************************************** + ! public subroutine computJacobSundials: compute the Jacobian matrix + ! ********************************************************************************************************** + subroutine computJacobSundials(& + ! input: model control + cj, & ! intent(in): this scalar changes whenever the step size or method order changes + dt, & ! intent(in): length of the time step (seconds) + nSnow, & ! intent(in): number of snow layers + nSoil, & ! intent(in): number of soil layers + nLayers, & ! intent(in): total number of layers + computeVegFlux, & ! intent(in): flag to indicate if we need to compute fluxes over vegetation + computeBaseflow, & ! intent(in): flag to indicate if we need to compute baseflow + ixMatrix, & ! intent(in): form of the Jacobian matrix + specificStorage, & ! intent(in): specific storage coefficient (m-1) + theta_sat, & ! intent(in): soil porosity (-) + ixRichards, & ! intent(in): choice of option for Richards' equation + ! input: data structures + indx_data, & ! intent(in): index data + prog_data, & ! intent(in): model prognostic variables for a local HRU + diag_data, & ! intent(in): model diagnostic variables for a local HRU + deriv_data, & ! intent(in): derivatives in model fluxes w.r.t. relevant state variables + dBaseflow_dMatric, & ! intent(in): derivative in baseflow w.r.t. matric head (s-1) + ! input: state variables + mLayerTempPrime, & ! intent(in): vector of derivative value for layer temperature (K) + mLayerMatricHeadPrime, & ! intent(in) + mLayerMatricHeadLiqPrime, & ! intent(in) + mLayerVolFracWatPrime, & ! intent(in) + scalarCanopyTempPrime, & ! intent(in): derivative value for temperature of the vegetation canopy (K) + scalarCanopyWatPrime, & ! intent(in) + ! input-output: Jacobian and its diagonal + dMat, & ! intent(inout): diagonal of the Jacobian matrix + aJac, & ! intent(out): Jacobian matrix + ! output: error control + err,message) ! intent(out): error code and error message + ! ----------------------------------------------------------------------------------------------------------------- + implicit none + ! input: model control + real(rkind),intent(in) :: cj + real(rkind),intent(in) :: dt ! length of the time step (seconds) + integer(i4b),intent(in) :: nSnow ! number of snow layers + integer(i4b),intent(in) :: nSoil ! number of soil layers + integer(i4b),intent(in) :: nLayers ! total number of layers in the snow+soil domain + logical(lgt),intent(in) :: computeVegFlux ! flag to indicate if computing fluxes over vegetation + logical(lgt),intent(in) :: computeBaseflow ! flag to indicate if computing baseflow + integer(i4b),intent(in) :: ixMatrix ! form of the Jacobian matrix + real(rkind),intent(in) :: specificStorage ! specific storage coefficient (m-1) + real(rkind),intent(in) :: theta_sat(:) ! soil porosity (-) + integer(i4b),intent(in) :: ixRichards ! choice of option for Richards' equation + ! input: data structures + type(var_ilength),intent(in) :: indx_data ! indices defining model states and layers + type(var_dlength),intent(in) :: prog_data ! prognostic variables for a local HRU + type(var_dlength),intent(in) :: diag_data ! diagnostic variables for a local HRU + type(var_dlength),intent(in) :: deriv_data ! derivatives in model fluxes w.r.t. relevant state variables + real(rkind),intent(in) :: dBaseflow_dMatric(:,:) ! derivative in baseflow w.r.t. matric head (s-1) + ! input: state variables + real(rkind),intent(in) :: mLayerTempPrime(:) + real(rkind),intent(in) :: mLayerMatricHeadPrime(:) + real(rkind),intent(in) :: mLayerMatricHeadLiqPrime(:) + real(rkind),intent(in) :: mLayerVolFracWatPrime(:) + real(rkind),intent(in) :: scalarCanopyTempPrime ! derivative value for temperature of the vegetation canopy (K) + real(rkind),intent(in) :: scalarCanopyWatPrime + ! input-output: Jacobian and its diagonal + real(rkind),intent(inout) :: dMat(:) ! diagonal of the Jacobian matrix + real(rkind),intent(out) :: aJac(:,:) ! Jacobian matrix + ! output variables + integer(i4b),intent(out) :: err ! error code + character(*),intent(out) :: message ! error message + ! -------------------------------------------------------------- + ! * local variables + ! -------------------------------------------------------------- + ! indices of model state variables + integer(i4b) :: jState ! index of state within the state subset + integer(i4b) :: qState ! index of cross-derivative state variable for baseflow + integer(i4b) :: nrgState ! energy state variable + integer(i4b) :: watState ! hydrology state variable + integer(i4b) :: nState ! number of state variables + ! indices of model layers + integer(i4b) :: iLayer ! index of model layer + integer(i4b) :: jLayer ! index of model layer within the full state vector (hydrology) + integer(i4b) :: pLayer ! indices of soil layers (used for the baseflow derivatives) + ! conversion factors + real(rkind) :: convLiq2tot ! factor to convert liquid water derivative to total water derivative + ! -------------------------------------------------------------- + ! associate variables from data structures + associate(& + ! indices of model state variables + ixCasNrg => indx_data%var(iLookINDEX%ixCasNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy air space energy state variable + ixVegNrg => indx_data%var(iLookINDEX%ixVegNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy energy state variable + ixVegHyd => indx_data%var(iLookINDEX%ixVegHyd)%dat(1) ,& ! intent(in): [i4b] index of canopy hydrology state variable (mass) + ixTopNrg => indx_data%var(iLookINDEX%ixTopNrg)%dat(1) ,& ! intent(in): [i4b] index of upper-most energy state in the snow+soil subdomain + ixTopHyd => indx_data%var(iLookINDEX%ixTopHyd)%dat(1) ,& ! intent(in): [i4b] index of upper-most hydrology state in the snow+soil subdomain + ixAqWat => indx_data%var(iLookINDEX%ixAqWat)%dat(1) ,& ! intent(in): [i4b] index of water storage in the aquifer + ! vectors of indices for specfic state types within specific sub-domains IN THE FULL STATE VECTOR + ixNrgLayer => indx_data%var(iLookINDEX%ixNrgLayer)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in the snow+soil domain + ixHydLayer => indx_data%var(iLookINDEX%ixHydLayer)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for hydrology states in the snow+soil domain + ! vector of energy indices for the snow and soil domains + ! NOTE: states not in the subset are equal to integerMissing + ixSnowSoilNrg => indx_data%var(iLookINDEX%ixSnowSoilNrg)%dat ,& ! intent(in): [i4b(:)] index in the state subset for energy state variables in the snow+soil domain + ixSnowOnlyNrg => indx_data%var(iLookINDEX%ixSnowOnlyNrg)%dat ,& ! intent(in): [i4b(:)] index in the state subset for energy state variables in the snow domain + ixSoilOnlyNrg => indx_data%var(iLookINDEX%ixSoilOnlyNrg)%dat ,& ! intent(in): [i4b(:)] index in the state subset for energy state variables in the soil domain + ! vector of hydrology indices for the snow and soil domains + ! NOTE: states not in the subset are equal to integerMissing + ixSnowSoilHyd => indx_data%var(iLookINDEX%ixSnowSoilHyd)%dat ,& ! intent(in): [i4b(:)] index in the state subset for hydrology state variables in the snow+soil domain + ixSnowOnlyHyd => indx_data%var(iLookINDEX%ixSnowOnlyHyd)%dat ,& ! intent(in): [i4b(:)] index in the state subset for hydrology state variables in the snow domain + ixSoilOnlyHyd => indx_data%var(iLookINDEX%ixSoilOnlyHyd)%dat ,& ! intent(in): [i4b(:)] index in the state subset for hydrology state variables in the soil domain + ! number of state variables of a specific type + nSnowSoilNrg => indx_data%var(iLookINDEX%nSnowSoilNrg )%dat(1) ,& ! intent(in): [i4b] number of energy state variables in the snow+soil domain + nSnowOnlyNrg => indx_data%var(iLookINDEX%nSnowOnlyNrg )%dat(1) ,& ! intent(in): [i4b] number of energy state variables in the snow domain + nSoilOnlyNrg => indx_data%var(iLookINDEX%nSoilOnlyNrg )%dat(1) ,& ! intent(in): [i4b] number of energy state variables in the soil domain + nSnowSoilHyd => indx_data%var(iLookINDEX%nSnowSoilHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the snow+soil domain + nSnowOnlyHyd => indx_data%var(iLookINDEX%nSnowOnlyHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the snow domain + nSoilOnlyHyd => indx_data%var(iLookINDEX%nSoilOnlyHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the soil domain + ! type and index of model control volume + ixHydType => indx_data%var(iLookINDEX%ixHydType)%dat ,& ! intent(in): [i4b(:)] index of the type of hydrology states in snow+soil domain + ixDomainType => indx_data%var(iLookINDEX%ixDomainType)%dat ,& ! intent(in): [i4b(:)] indices defining the type of the domain (iname_veg, iname_snow, iname_soil) + ixControlVolume => indx_data%var(iLookINDEX%ixControlVolume)%dat ,& ! intent(in): [i4b(:)] index of the control volume for specific model domains + ! mapping between states and model layers + ixMapSubset2Full => indx_data%var(iLookINDEX%ixMapSubset2Full)%dat ,& ! intent(in): [i4b(:)] list of indices in the full state vector that are in the state subset + ixMapFull2Subset => indx_data%var(iLookINDEX%ixMapFull2Subset)%dat ,& ! intent(in): [i4b(:)] list of indices in the state subset in each element of the full state vector + ! derivatives in net vegetation energy fluxes w.r.t. relevant state variables + dCanairNetFlux_dCanairTemp => deriv_data%var(iLookDERIV%dCanairNetFlux_dCanairTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy air space flux w.r.t. canopy air temperature + dCanairNetFlux_dCanopyTemp => deriv_data%var(iLookDERIV%dCanairNetFlux_dCanopyTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy air space flux w.r.t. canopy temperature + dCanairNetFlux_dGroundTemp => deriv_data%var(iLookDERIV%dCanairNetFlux_dGroundTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy air space flux w.r.t. ground temperature + dCanopyNetFlux_dCanairTemp => deriv_data%var(iLookDERIV%dCanopyNetFlux_dCanairTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy flux w.r.t. canopy air temperature + dCanopyNetFlux_dCanopyTemp => deriv_data%var(iLookDERIV%dCanopyNetFlux_dCanopyTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy flux w.r.t. canopy temperature + dCanopyNetFlux_dGroundTemp => deriv_data%var(iLookDERIV%dCanopyNetFlux_dGroundTemp )%dat(1) ,& ! intent(in): [dp] derivative in net canopy flux w.r.t. ground temperature + dCanopyNetFlux_dCanWat => deriv_data%var(iLookDERIV%dCanopyNetFlux_dCanWat )%dat(1) ,& ! intent(in): [dp] derivative in net canopy fluxes w.r.t. canopy total water content + dGroundNetFlux_dCanairTemp => deriv_data%var(iLookDERIV%dGroundNetFlux_dCanairTemp )%dat(1) ,& ! intent(in): [dp] derivative in net ground flux w.r.t. canopy air temperature + dGroundNetFlux_dCanopyTemp => deriv_data%var(iLookDERIV%dGroundNetFlux_dCanopyTemp )%dat(1) ,& ! intent(in): [dp] derivative in net ground flux w.r.t. canopy temperature + dGroundNetFlux_dCanWat => deriv_data%var(iLookDERIV%dGroundNetFlux_dCanWat )%dat(1) ,& ! intent(in): [dp] derivative in net ground fluxes w.r.t. canopy total water content + ! derivatives in evaporative fluxes w.r.t. relevant state variables + dCanopyEvaporation_dTCanair => deriv_data%var(iLookDERIV%dCanopyEvaporation_dTCanair )%dat(1) ,& ! intent(in): [dp] derivative in canopy evaporation w.r.t. canopy air temperature + dCanopyEvaporation_dTCanopy => deriv_data%var(iLookDERIV%dCanopyEvaporation_dTCanopy )%dat(1) ,& ! intent(in): [dp] derivative in canopy evaporation w.r.t. canopy temperature + dCanopyEvaporation_dTGround => deriv_data%var(iLookDERIV%dCanopyEvaporation_dTGround )%dat(1) ,& ! intent(in): [dp] derivative in canopy evaporation w.r.t. ground temperature + dCanopyEvaporation_dCanWat => deriv_data%var(iLookDERIV%dCanopyEvaporation_dCanWat )%dat(1) ,& ! intent(in): [dp] derivative in canopy evaporation w.r.t. canopy total water content + dGroundEvaporation_dTCanair => deriv_data%var(iLookDERIV%dGroundEvaporation_dTCanair )%dat(1) ,& ! intent(in): [dp] derivative in ground evaporation w.r.t. canopy air temperature + dGroundEvaporation_dTCanopy => deriv_data%var(iLookDERIV%dGroundEvaporation_dTCanopy )%dat(1) ,& ! intent(in): [dp] derivative in ground evaporation w.r.t. canopy temperature + dGroundEvaporation_dTGround => deriv_data%var(iLookDERIV%dGroundEvaporation_dTGround )%dat(1) ,& ! intent(in): [dp] derivative in ground evaporation w.r.t. ground temperature + dGroundEvaporation_dCanWat => deriv_data%var(iLookDERIV%dGroundEvaporation_dCanWat )%dat(1) ,& ! intent(in): [dp] derivative in ground evaporation w.r.t. canopy total water content + ! derivatives in canopy water w.r.t canopy temperature + dCanLiq_dTcanopy => deriv_data%var(iLookDERIV%dCanLiq_dTcanopy )%dat(1) ,& ! intent(in): [dp] derivative in canopy liquid storage w.r.t. temperature + dTheta_dTkCanopy => deriv_data%var(iLookDERIV%dTheta_dTkCanopy )%dat(1) ,& ! intent(in): [dp] derivative in volumetric liquid water content w.r.t. temperature + d2Theta_dTkCanopy2 => deriv_data%var(iLookDERIV%d2Theta_dTkCanopy2 )%dat(1) ,& ! intent(in): [dp] second derivative of volumetric liquid water content w.r.t. temperature + dFracLiqVeg_dTkCanopy => deriv_data%var(iLookDERIV%dFracLiqVeg_dTkCanopy )%dat(1) ,& ! intent(in): [dp] derivative in fraction of (throughfall + drainage) w.r.t. temperature + ! derivatives in canopy liquid fluxes w.r.t. canopy water + scalarCanopyLiqDeriv => deriv_data%var(iLookDERIV%scalarCanopyLiqDeriv )%dat(1) ,& ! intent(in): [dp] derivative in (throughfall + drainage) w.r.t. canopy liquid water + ! derivatives in energy fluxes at the interface of snow+soil layers w.r.t. temperature in layers above and below + dNrgFlux_dTempAbove => deriv_data%var(iLookDERIV%dNrgFlux_dTempAbove )%dat ,& ! intent(in): [dp(:)] derivatives in the flux w.r.t. temperature in the layer above + dNrgFlux_dTempBelow => deriv_data%var(iLookDERIV%dNrgFlux_dTempBelow )%dat ,& ! intent(in): [dp(:)] derivatives in the flux w.r.t. temperature in the layer below + ! derivatives in energy fluxes at the interface of snow+soil layers w.r.t. water state in layers above and below + dNrgFlux_dWatAbove => deriv_data%var(iLookDERIV%dNrgFlux_dWatAbove )%dat ,& ! intent(in): [dp(:)] derivatives in the flux w.r.t. water state in the layer above + dNrgFlux_dWatBelow => deriv_data%var(iLookDERIV%dNrgFlux_dWatBelow )%dat ,& ! intent(in): [dp(:)] derivatives in the flux w.r.t. water state in the layer below + ! derivatives in soil transpiration w.r.t. canopy state variables + mLayerdTrans_dTCanair => deriv_data%var(iLookDERIV%mLayerdTrans_dTCanair )%dat ,& ! intent(in): derivatives in the soil layer transpiration flux w.r.t. canopy air temperature + mLayerdTrans_dTCanopy => deriv_data%var(iLookDERIV%mLayerdTrans_dTCanopy )%dat ,& ! intent(in): derivatives in the soil layer transpiration flux w.r.t. canopy temperature + mLayerdTrans_dTGround => deriv_data%var(iLookDERIV%mLayerdTrans_dTGround )%dat ,& ! intent(in): derivatives in the soil layer transpiration flux w.r.t. ground temperature + mLayerdTrans_dCanWat => deriv_data%var(iLookDERIV%mLayerdTrans_dCanWat )%dat ,& ! intent(in): derivatives in the soil layer transpiration flux w.r.t. canopy total water + ! derivatives in aquifer transpiration w.r.t. canopy state variables + dAquiferTrans_dTCanair => deriv_data%var(iLookDERIV%dAquiferTrans_dTCanair )%dat(1) ,& ! intent(in): derivatives in the aquifer transpiration flux w.r.t. canopy air temperature + dAquiferTrans_dTCanopy => deriv_data%var(iLookDERIV%dAquiferTrans_dTCanopy )%dat(1) ,& ! intent(in): derivatives in the aquifer transpiration flux w.r.t. canopy temperature + dAquiferTrans_dTGround => deriv_data%var(iLookDERIV%dAquiferTrans_dTGround )%dat(1) ,& ! intent(in): derivatives in the aquifer transpiration flux w.r.t. ground temperature + dAquiferTrans_dCanWat => deriv_data%var(iLookDERIV%dAquiferTrans_dCanWat )%dat(1) ,& ! intent(in): derivatives in the aquifer transpiration flux w.r.t. canopy total water + ! derivative in liquid water fluxes at the interface of snow layers w.r.t. volumetric liquid water content in the layer above + iLayerLiqFluxSnowDeriv => deriv_data%var(iLookDERIV%iLayerLiqFluxSnowDeriv )%dat ,& ! intent(in): [dp(:)] derivative in vertical liquid water flux at layer interfaces + ! derivative in liquid water fluxes for the soil domain w.r.t hydrology state variables + dVolTot_dPsi0 => deriv_data%var(iLookDERIV%dVolTot_dPsi0 )%dat ,& ! intent(in): [dp(:)] derivative in total water content w.r.t. total water matric potential + d2VolTot_d2Psi0 => deriv_data%var(iLookDERIV%d2VolTot_d2Psi0 )%dat ,& ! intent(in): [dp(:)] second derivative in total water content w.r.t. total water matric potential + dCompress_dPsi => deriv_data%var(iLookDERIV%dCompress_dPsi )%dat ,& ! intent(in): [dp(:)] derivative in compressibility w.r.t matric head + dq_dHydStateAbove => deriv_data%var(iLookDERIV%dq_dHydStateAbove )%dat ,& ! intent(in): [dp(:)] change in flux at layer interfaces w.r.t. states in the layer above + dq_dHydStateBelow => deriv_data%var(iLookDERIV%dq_dHydStateBelow )%dat ,& ! intent(in): [dp(:)] change in flux at layer interfaces w.r.t. states in the layer below + dq_dHydStateLayerSurfVec => deriv_data%var(iLookDERIV%dq_dHydStateLayerSurfVec )%dat ,& ! intent(in): [dp(:)] change in the flux in soil surface interface w.r.t. state variables in layers + ! derivative in baseflow flux w.r.t. aquifer storage + dBaseflow_dAquifer => deriv_data%var(iLookDERIV%dBaseflow_dAquifer )%dat(1) ,& ! intent(in): [dp(:)] derivative in baseflow flux w.r.t. aquifer storage (s-1) + ! derivative in liquid water fluxes for the soil domain w.r.t energy state variables + dq_dNrgStateAbove => deriv_data%var(iLookDERIV%dq_dNrgStateAbove )%dat ,& ! intent(in): [dp(:)] change in flux at layer interfaces w.r.t. states in the layer above + dq_dNrgStateBelow => deriv_data%var(iLookDERIV%dq_dNrgStateBelow )%dat ,& ! intent(in): [dp(:)] change in flux at layer interfaces w.r.t. states in the layer below + dq_dNrgStateLayerSurfVec => deriv_data%var(iLookDERIV%dq_dNrgStateLayerSurfVec )%dat ,& ! intent(in): [dp(:)] change in the flux in soil surface interface w.r.t. state variables in layers + ! derivative in liquid water fluxes for the soil and snow domain w.r.t temperature + dFracLiqSnow_dTk => deriv_data%var(iLookDERIV%dFracLiqSnow_dTk )%dat ,& ! intent(in): [dp(:)] derivative in fraction of liquid snow w.r.t. temperature + mLayerdTheta_dTk => deriv_data%var(iLookDERIV%mLayerdTheta_dTk )%dat ,& ! intent(in): [dp(:)] derivative in volumetric liquid water content w.r.t. temperature + mLayerd2Theta_dTk2 => deriv_data%var(iLookDERIV%mLayerd2Theta_dTk2 )%dat ,& ! intent(in): [dp(:)] second derivative of volumetric liquid water content w.r.t. temperature + ! derivate in bulk heat capacity w.r.t. relevant state variables + dVolHtCapBulk_dPsi0 => deriv_data%var(iLookDERIV%dVolHtCapBulk_dPsi0 )%dat ,& ! intent(in): [dp(:)] derivative in bulk heat capacity w.r.t. matric potential + dVolHtCapBulk_dTheta => deriv_data%var(iLookDERIV%dVolHtCapBulk_dTheta )%dat ,& ! intent(in): [dp(:)] derivative in bulk heat capacity w.r.t. volumetric water content + dVolHtCapBulk_dCanWat => deriv_data%var(iLookDERIV%dVolHtCapBulk_dCanWat )%dat(1) ,& ! intent(in): [dp] derivative in bulk heat capacity w.r.t. volumetric water content + dVolHtCapBulk_dTk => deriv_data%var(iLookDERIV%dVolHtCapBulk_dTk )%dat ,& ! intent(in): [dp(:)] derivative in bulk heat capacity w.r.t. temperature + dVolHtCapBulk_dTkCanopy => deriv_data%var(iLookDERIV%dVolHtCapBulk_dTkCanopy )%dat(1) ,& ! intent(in): [dp] derivative in bulk heat capacity w.r.t. temperature + ! diagnostic variables + scalarFracLiqVeg => diag_data%var(iLookDIAG%scalarFracLiqVeg )%dat(1) ,& ! intent(in): [dp] fraction of liquid water on vegetation (-) + scalarBulkVolHeatCapVeg => diag_data%var(iLookDIAG%scalarBulkVolHeatCapVeg )%dat(1) ,& ! intent(in): [dp] bulk volumetric heat capacity of vegetation (J m-3 K-1) + mLayerFracLiqSnow => diag_data%var(iLookDIAG%mLayerFracLiqSnow )%dat ,& ! intent(in): [dp(:)] fraction of liquid water in each snow layer (-) + mLayerVolHtCapBulk => diag_data%var(iLookDIAG%mLayerVolHtCapBulk )%dat ,& ! intent(in): [dp(:)] bulk volumetric heat capacity in each snow and soil layer (J m-3 K-1) + scalarSoilControl => diag_data%var(iLookDIAG%scalarSoilControl )%dat(1) ,& ! intent(in): [dp] soil control on infiltration, zero or one + ! canopy and layer depth + canopyDepth => diag_data%var(iLookDIAG%scalarCanopyDepth )%dat(1) ,& ! intent(in): [dp ] canopy depth (m) + mLayerDepth => prog_data%var(iLookPROG%mLayerDepth )%dat ,& ! intent(in): [dp(:)] depth of each layer in the snow-soil sub-domain (m) + layerType => indx_data%var(iLookINDEX%layerType )%dat & ! intent(in): [i4b(:)] named variables defining the type of layer in snow+soil domain + ) ! making association with data in structures + ! -------------------------------------------------------------- + ! initialize error control + err=0; message='computJacobSundials/' + + ! ********************************************************************************************************************************************************* + ! ********************************************************************************************************************************************************* + ! * PART 0: PRELIMINARIES (INITIALIZE JACOBIAN AND COMPUTE TIME-VARIABLE DIAGONAL TERMS) + ! ********************************************************************************************************************************************************* + ! ********************************************************************************************************************************************************* + + ! get the number of state variables + nState = size(dMat) + + ! initialize the Jacobian + ! NOTE: this needs to be done every time, since Jacobian matrix is modified in the solver + aJac(:,:) = 0._rkind ! analytical Jacobian matrix + + ! compute terms in the Jacobian for vegetation (excluding fluxes) + ! NOTE: energy for vegetation is computed *within* the iteration loop as it includes phase change + if(ixVegNrg/=integerMissing)then + dMat(ixVegNrg) = ( scalarBulkVolHeatCapVeg + LH_fus*iden_water*dTheta_dTkCanopy ) * cj & + + dVolHtCapBulk_dTkCanopy * scalarCanopyTempPrime & + + LH_fus*iden_water * scalarCanopyTempPrime * d2Theta_dTkCanopy2 & + + LH_fus * dFracLiqVeg_dTkCanopy * scalarCanopyWatPrime / canopyDepth + endif + + ! compute additional terms for the Jacobian for the snow-soil domain (excluding fluxes) + ! NOTE: energy for snow+soil is computed *within* the iteration loop as it includes phase change + do iLayer=1,nLayers + if(ixSnowSoilNrg(iLayer)/=integerMissing)then + dMat(ixSnowSoilNrg(iLayer)) = ( mLayerVolHtCapBulk(iLayer) + LH_fus*iden_water*mLayerdTheta_dTk(iLayer) ) * cj & + + dVolHtCapBulk_dTk(iLayer) * mLayerTempPrime(iLayer) & + + LH_fus*iden_water * mLayerTempPrime(iLayer) * mLayerd2Theta_dTk2(iLayer) & + + LH_fus*iden_water * dFracLiqSnow_dTk(iLayer) * mLayerVolFracWatPrime(iLayer) + endif + end do + + ! compute additional terms for the Jacobian for the soil domain (excluding fluxes) + do iLayer=1,nSoil + if(ixSoilOnlyHyd(iLayer)/=integerMissing)then + dMat(ixSoilOnlyHyd(iLayer)) = ( dVolTot_dPsi0(iLayer) + dCompress_dPsi(iLayer) ) * cj + d2VolTot_d2Psi0(iLayer) * mLayerMatricHeadPrime(iLayer) + + if(ixRichards==mixdform)then + dMat(ixSoilOnlyHyd(iLayer)) = dMat(ixSoilOnlyHyd(iLayer)) + specificStorage * dVolTot_dPsi0(iLayer) * mLayerMatricHeadPrime(iLayer) / theta_sat(iLayer) + endif + + endif + end do + + ! define the form of the matrix + select case(ixMatrix) + + ! ********************************************************************************************************************************************************* + ! ********************************************************************************************************************************************************* + ! * PART 1: BAND MATRIX + ! ********************************************************************************************************************************************************* + ! ********************************************************************************************************************************************************* + case(ixBandMatrix) ! ixBandMatrix ixFullMatrix + print *, 'banded jacobian matrix needs to be implemented' + stop 1 + + ! ********************************************************************************************************************************************************* + ! ********************************************************************************************************************************************************* + ! * PART 2: FULL MATRIX + ! ********************************************************************************************************************************************************* + ! ********************************************************************************************************************************************************* + case(ixFullMatrix) ! ixFullMatrix ixBandMatrix + + ! check + if(size(aJac,1)/=size(dMat) .or. size(aJac,2)/=size(dMat))then + message=trim(message)//'unexpected shape of the Jacobian matrix: expect aJac(nState,nState)' + err=20; return + endif + + ! ----- + ! * energy and liquid fluxes over vegetation... + ! --------------------------------------------- + if(computeVegFlux)then ! (derivatives only defined when vegetation protrudes over the surface) + + ! * energy fluxes with the canopy water + if(ixVegHyd/=integerMissing)then + + ! * cross-derivative terms w.r.t. system temperatures (kg m-2 K-1) + if(ixCasNrg/=integerMissing) aJac(ixVegHyd,ixCasNrg) = -dCanopyEvaporation_dTCanair*dt + ! dt*scalarCanopyLiqDeriv*dCanLiq_dTcanopy is the derivative in throughfall and canopy drainage with canopy temperature + if(ixVegNrg/=integerMissing) aJac(ixVegHyd,ixVegNrg) = -dCanopyEvaporation_dTCanopy*dt + dt*scalarCanopyLiqDeriv*dCanLiq_dTcanopy + ! * liquid water fluxes for vegetation canopy (-), dt*scalarFracLiqVeg*scalarCanopyLiqDeriv is the derivative in throughfall and canopy drainage with canopy water + aJac(ixVegHyd,ixVegHyd) = -scalarFracLiqVeg*(dCanopyEvaporation_dCanWat - scalarCanopyLiqDeriv)*dt + 1._rkind * cj + if(ixTopNrg/=integerMissing) aJac(ixVegHyd,ixTopNrg) = -dCanopyEvaporation_dTGround*dt + + ! * cross-derivative terms w.r.t. canopy water (kg-1 m2) + if(ixTopHyd/=integerMissing) aJac(ixTopHyd,ixVegHyd) = (dt/mLayerDepth(1))*(-scalarSoilControl*scalarFracLiqVeg*scalarCanopyLiqDeriv)/iden_water + + ! * cross-derivative terms w.r.t. canopy liquid water (J m-1 kg-1) + ! NOTE: dIce/dLiq = (1 - scalarFracLiqVeg); dIce*LH_fus/canopyDepth = J m-3; dLiq = kg m-2 + if(ixVegNrg/=integerMissing) aJac(ixVegNrg,ixVegHyd) = (-1._rkind + scalarFracLiqVeg)*LH_fus/canopyDepth * cj & + + dVolHtCapBulk_dCanWat * scalarCanopyTempPrime - (dt/canopyDepth) * dCanopyNetFlux_dCanWat & + + LH_fus * scalarCanopyTempPrime * dFracLiqVeg_dTkCanopy / canopyDepth ! dF/dLiq + if(ixTopNrg/=integerMissing) aJac(ixTopNrg,ixVegHyd) = (dt/mLayerDepth(1))*(-dGroundNetFlux_dCanWat) + endif + + ! * -derivative terms w.r.t. canopy temperature (K-1) + if(ixVegNrg/=integerMissing)then + if(ixTopHyd/=integerMissing) aJac(ixTopHyd,ixVegNrg) = (dt/mLayerDepth(1))*(-scalarSoilControl*scalarCanopyLiqDeriv*dCanLiq_dTcanopy)/iden_water + endif + + ! * energy fluxes with the canopy air space (J m-3 K-1) + if(ixCasNrg/=integerMissing)then + aJac(ixCasNrg,ixCasNrg) = (dt/canopyDepth)*(-dCanairNetFlux_dCanairTemp) + dMat(ixCasNrg) * cj + if(ixVegNrg/=integerMissing) aJac(ixCasNrg,ixVegNrg) = (dt/canopyDepth)*(-dCanairNetFlux_dCanopyTemp) + if(ixTopNrg/=integerMissing) aJac(ixCasNrg,ixTopNrg) = (dt/canopyDepth)*(-dCanairNetFlux_dGroundTemp) + endif + + ! * energy fluxes with the vegetation canopy (J m-3 K-1) + if(ixVegNrg/=integerMissing)then + if(ixCasNrg/=integerMissing) aJac(ixVegNrg,ixCasNrg) = (dt/canopyDepth)*(-dCanopyNetFlux_dCanairTemp) + aJac(ixVegNrg,ixVegNrg) = (dt/canopyDepth)*(-dCanopyNetFlux_dCanopyTemp) + dMat(ixVegNrg) + if(ixTopNrg/=integerMissing) aJac(ixVegNrg,ixTopNrg) = (dt/canopyDepth)*(-dCanopyNetFlux_dGroundTemp) + endif + + ! * energy fluxes with the surface (J m-3 K-1) + if(ixTopNrg/=integerMissing)then + if(ixCasNrg/=integerMissing) aJac(ixTopNrg,ixCasNrg) = (dt/mLayerDepth(1))*(-dGroundNetFlux_dCanairTemp) + if(ixVegNrg/=integerMissing) aJac(ixTopNrg,ixVegNrg) = (dt/mLayerDepth(1))*(-dGroundNetFlux_dCanopyTemp) + endif + + endif ! if there is a need to compute energy fluxes within vegetation + + ! ----- + ! * energy fluxes for the snow+soil domain... + ! ------------------------------------------- + if(nSnowSoilNrg>0)then + do iLayer=1,nLayers ! loop through all layers in the snow+soil domain + + ! check if the state is in the subset + if(ixSnowSoilNrg(iLayer)==integerMissing) cycle + + ! - define index within the state subset and the full state vector + jState = ixSnowSoilNrg(iLayer) ! index within the state subset + + ! - diagonal elements + aJac(jState,jState) = (dt/mLayerDepth(iLayer))*(-dNrgFlux_dTempBelow(iLayer-1) + dNrgFlux_dTempAbove(iLayer)) + dMat(jState) + + ! - lower-diagonal elements + if(iLayer>1)then + if(ixSnowSoilNrg(iLayer-1)/=integerMissing) aJac(ixSnowSoilNrg(iLayer-1),jState) = (dt/mLayerDepth(iLayer-1))*( dNrgFlux_dTempBelow(iLayer-1) ) + endif + + ! - upper diagonal elements + if(iLayer<nLayers)then + if(ixSnowSoilNrg(iLayer+1)/=integerMissing) aJac(ixSnowSoilNrg(iLayer+1),jState) = (dt/mLayerDepth(iLayer+1))*(-dNrgFlux_dTempAbove(iLayer ) ) + endif + + end do ! (looping through energy states in the snow+soil domain) + endif ! (if the subset includes energy state variables in the snow+soil domain) + + ! ----- + ! * liquid water fluxes for the snow domain... + ! -------------------------------------------- + if(nSnowOnlyHyd>0)then + do iLayer=1,nSnow ! loop through layers in the snow domain + + ! - check that the snow layer is desired + if(ixSnowOnlyHyd(iLayer)==integerMissing) cycle + + ! - define state indices for the current layer + watState = ixSnowOnlyHyd(iLayer) ! hydrology state index within the state subset + + ! compute factor to convert liquid water derivative to total water derivative + select case( ixHydType(iLayer) ) + case(iname_watLayer); convLiq2tot = mLayerFracLiqSnow(iLayer) + case default; convLiq2tot = 1._rkind + end select + + ! - diagonal elements + aJac(watState,watState) = (dt/mLayerDepth(iLayer))*iLayerLiqFluxSnowDeriv(iLayer)*convLiq2tot + dMat(watState) * cj + + ! - lower-diagonal elements + if(iLayer>1)then + if(ixSnowOnlyHyd(iLayer-1)/=integerMissing) aJac(ixSnowOnlyHyd(iLayer-1),watState) = 0._rkind ! sub-diagonal: no dependence on other layers + endif + + ! - upper diagonal elements + if(iLayer<nSnow)then + if(ixSnowOnlyHyd(iLayer+1)/=integerMissing) aJac(ixSnowOnlyHyd(iLayer+1),watState) = -(dt/mLayerDepth(iLayer+1))*iLayerLiqFluxSnowDeriv(iLayer)*convLiq2tot ! dVol(below)/dLiq(above) -- (-) + endif + + end do ! (looping through liquid water states in the snow domain) + endif ! (if the subset includes hydrology state variables in the snow domain) + + ! ----- + ! * cross derivatives in the snow domain... + ! ---------------------------------------- + if(nSnowOnlyHyd>0 .and. nSnowOnlyNrg>0)then + do iLayer=1,nSnow ! loop through layers in the snow domain + + ! - check that the snow layer is desired + if(ixSnowOnlyNrg(iLayer)==integerMissing) cycle + + ! (define the energy state) + nrgState = ixSnowOnlyNrg(iLayer) ! index within the full state vector + + ! - define state indices for the current layer + watState = ixSnowOnlyHyd(iLayer) ! hydrology state index within the state subset + + if(watstate/=integerMissing)then ! (energy state for the current layer is within the state subset) + + ! - include derivatives of energy fluxes w.r.t water fluxes for current layer + aJac(nrgState,watState) = (-1._rkind + mLayerFracLiqSnow(iLayer))*LH_fus*iden_water * cj & + + dVolHtCapBulk_dTheta(iLayer) * mLayerTempPrime(iLayer) & + + (dt/mLayerDepth(iLayer))*(-dNrgFlux_dWatBelow(iLayer-1) + dNrgFlux_dWatAbove(iLayer)) & + + LH_fus*iden_water * mLayerTempPrime(iLayer) * dFracLiqSnow_dTk(iLayer) ! (dF/dLiq) + + ! - include derivatives of water fluxes w.r.t energy fluxes for current layer + aJac(watState,nrgState) = (dt/mLayerDepth(iLayer))*iLayerLiqFluxSnowDeriv(iLayer)*mLayerdTheta_dTk(iLayer) ! (dVol/dT) + + ! (cross-derivative terms for the layer below) + if(iLayer<nSnow)then + aJac(ixSnowOnlyHyd(iLayer+1),nrgState) = -(dt/mLayerDepth(iLayer+1))*iLayerLiqFluxSnowDeriv(iLayer)*mLayerdTheta_dTk(iLayer) ! dVol(below)/dT(above) -- K-1 + endif ! (if there is a water state in the layer below the current layer in the given state subset) + + ! - include derivatives of heat capacity w.r.t water fluxes for surrounding layers starting with layer above + if(iLayer>1)then + if(ixSnowOnlyNrg(iLayer-1)/=integerMissing) aJac(ixSnowOnlyNrg(iLayer-1),watState) = (dt/mLayerDepth(iLayer-1))*( dNrgFlux_dWatBelow(iLayer-1) ) + endif + + ! (cross-derivative terms for the layer below) + if(iLayer<nSnow)then + if(ixSnowOnlyNrg(iLayer+1)/=integerMissing) aJac(ixSnowOnlyNrg(iLayer+1),watState) = (dt/mLayerDepth(iLayer+1))*(-dNrgFlux_dWatAbove(iLayer ) ) + elseif(iLayer==nSnow .and. nSoilOnlyNrg>0)then !bottom snow layer and there is soil below + if(ixSoilOnlyNrg(1)/=integerMissing) aJac(ixSoilOnlyNrg(1),watState) = (dt/mLayerDepth(nSnow+1))*(-dNrgFlux_dWatAbove(nSnow) ) + endif + + endif ! (if the energy state for the current layer is within the state subset) + + end do ! (looping through snow layers) + endif ! (if there are state variables for both water and energy in the snow domain) + + ! ----- + ! * liquid water fluxes for the soil domain... + ! -------------------------------------------- + if(nSoilOnlyHyd>0)then + do iLayer=1,nSoil + + ! - check that the soil layer is desired + if(ixSoilOnlyHyd(iLayer)==integerMissing) cycle + + ! - define state indices + watState = ixSoilOnlyHyd(iLayer) ! hydrology state index within the state subset + + ! - define indices of the soil layers + jLayer = iLayer+nSnow ! index of layer in the snow+soil vector + + ! - compute the diagonal elements + ! all terms *excluding* baseflow + aJac(watState,watState) = (dt/mLayerDepth(jLayer))*(-dq_dHydStateBelow(iLayer-1) + dq_dHydStateAbove(iLayer)) + dMat(watState) + + ! - compute the lower-diagonal elements + if(iLayer>1)then + if(ixSoilOnlyHyd(iLayer-1)/=integerMissing) aJac(ixSoilOnlyHyd(iLayer-1),watState) = (dt/mLayerDepth(jLayer-1))*( dq_dHydStateBelow(iLayer-1)) + endif + + ! - compute the upper-diagonal elements + if(iLayer<nSoil)then + if(ixSoilOnlyHyd(iLayer+1)/=integerMissing) aJac(ixSoilOnlyHyd(iLayer+1),watState) = (dt/mLayerDepth(jLayer+1))*(-dq_dHydStateAbove(iLayer)) + endif + + ! - include terms for baseflow + if(computeBaseflow .and. nSoilOnlyHyd==nSoil)then + do pLayer=1,nSoil + qState = ixSoilOnlyHyd(pLayer) ! hydrology state index within the state subset + aJac(watState,qState) = (dt/mLayerDepth(jLayer))*dBaseflow_dMatric(iLayer,pLayer) + aJac(watState,qState) + end do + endif + + ! - include terms for surface infiltration below surface + if(ixSoilOnlyHyd(1)/=integerMissing) aJac(ixSoilOnlyHyd(1),watState) = -(dt/mLayerDepth(1+nSnow))*dq_dHydStateLayerSurfVec(iLayer) + aJac(ixSoilOnlyHyd(1),watState) + + end do ! (looping through hydrology states in the soil domain) + + ! - include terms for surface infiltration above surface + if(nSnowOnlyHyd>0 .and. ixSnowOnlyHyd(nSnow)/=integerMissing)then + if(ixSoilOnlyHyd(1)/=integerMissing) aJac(ixSoilOnlyHyd(1),ixSnowOnlyHyd(nSnow)) = -(dt/mLayerDepth(1+nSnow))*dq_dHydStateLayerSurfVec(0) + elseif(computeVegFlux .and. ixVegHyd/=integerMissing)then !ixTopHyd = ixSoilOnlyHyd(1) + if(ixTopHyd/=integerMissing) aJac(ixTopHyd,ixVegHyd) = -(dt/mLayerDepth(1+nSnow))*dq_dHydStateLayerSurfVec(0) + aJac(ixTopHyd,ixVegHyd) + endif + + endif ! (if the subset includes hydrology state variables in the soil domain) + + ! ----- + ! * liquid water fluxes for the aquifer... + ! ---------------------------------------- + if(ixAqWat/=integerMissing) then + aJac(ixAqWat,ixAqWat) = -dBaseflow_dAquifer*dt + dMat(ixAqWat) * cj + aJac(ixAqWat,ixSoilOnlyNrg(nSoil)) = -dq_dNrgStateAbove(nSoil)*dt ! dAquiferRecharge_dTk = d_iLayerLiqFluxSoil(nSoil)_dTk + aJac(ixAqWat,ixSoilOnlyHyd(nSoil)) = -dq_dHydStateAbove(nSoil)*dt ! dAquiferRecharge_dWat = d_iLayerLiqFluxSoil(nSoil)_dWat + ! - include derivatives of energy and water w.r.t soil transpiration (dependent on canopy transpiration) + if(computeVegFlux)then + aJac(ixAqWat,ixCasNrg) = -dAquiferTrans_dTCanair*dt ! dVol/dT (K-1) + aJac(ixAqWat,ixVegNrg) = -dAquiferTrans_dTCanopy*dt ! dVol/dT (K-1) + aJac(ixAqWat,ixVegHyd) = -dAquiferTrans_dCanWat*dt ! dVol/dLiq (kg m-2)-1 + aJac(ixAqWat,ixTopNrg) = -dAquiferTrans_dTGround*dt ! dVol/dT (K-1) + endif + endif + + ! ----- + ! * cross derivatives in the soil domain... + ! ---------------------------------------- + if(nSoilOnlyHyd>0 .and. nSoilOnlyNrg>0)then + do iLayer=1,nSoilOnlyNrg + + ! - check that the soil layer is desired + if(ixSoilOnlyNrg(iLayer)==integerMissing) cycle + + ! - define indices of the soil layers + jLayer = iLayer+nSnow ! index of layer in the snow+soil vector + + ! - define the energy state variable + nrgState = ixNrgLayer(jLayer) ! index within the full state vector + + ! - define index of hydrology state variable within the state subset + watState = ixSoilOnlyHyd(iLayer) + + ! only compute derivatives if the water state for the current layer is within the state subset + if(watstate/=integerMissing)then + + ! - include derivatives in liquid water fluxes w.r.t. temperature for current layer + aJac(watState,nrgState) = (dt/mLayerDepth(jLayer))*(-dq_dNrgStateBelow(iLayer-1) + dq_dNrgStateAbove(iLayer)) ! dVol/dT (K-1) -- flux depends on ice impedance + + ! - compute lower diagonal elements + if(iLayer>1)then + if(ixSoilOnlyHyd(iLayer-1)/=integerMissing) aJac(ixSoilOnlyHyd(iLayer-1),nrgState) = (dt/mLayerDepth(jLayer-1))*( dq_dNrgStateBelow(iLayer-1)) ! K-1 + endif + + ! compute upper-diagonal elements + if(iLayer<nSoil)then + if(ixSoilOnlyHyd(iLayer+1)/=integerMissing) aJac(ixSoilOnlyHyd(iLayer+1),nrgState) = (dt/mLayerDepth(jLayer+1))*(-dq_dNrgStateAbove(iLayer)) ! K-1 + endif + + ! - include derivatives of energy w.r.t. ground evaporation + if(nSnow==0 .and. iLayer==1)then ! upper-most soil layer + if(computeVegFlux)then + aJac(watState,ixVegHyd) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dCanWat/iden_water) ! dVol/dLiq (kg m-2)-1 + aJac(watState,ixCasNrg) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dTCanair/iden_water) ! dVol/dT (K-1) + aJac(watState,ixVegNrg) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dTCanopy/iden_water) ! dVol/dT (K-1) + endif + aJac(watState,ixTopNrg) = (dt/mLayerDepth(jLayer))*(-dGroundEvaporation_dTGround/iden_water) + aJac(watState,ixTopNrg) ! dVol/dT (K-1) + endif + + ! - include derivatives of energy and water w.r.t soil transpiration (dependent on canopy transpiration) + if(computeVegFlux)then + aJac(watState,ixCasNrg) = (dt/mLayerDepth(jLayer))*(-mLayerdTrans_dTCanair(iLayer)) + aJac(watState,ixCasNrg) ! dVol/dT (K-1) + aJac(watState,ixVegNrg) = (dt/mLayerDepth(jLayer))*(-mLayerdTrans_dTCanopy(iLayer)) + aJac(watState,ixVegNrg) ! dVol/dT (K-1) + aJac(watState,ixVegHyd) = (dt/mLayerDepth(jLayer))*(-mLayerdTrans_dCanWat(iLayer)) + aJac(watState,ixVegHyd) ! dVol/dLiq (kg m-2)-1 + aJac(watState,ixTopNrg) = (dt/mLayerDepth(jLayer))*(-mLayerdTrans_dTGround(iLayer)) + aJac(watState,ixTopNrg) ! dVol/dT (K-1) + endif + + ! - include derivatives in energy fluxes w.r.t. with respect to water for current layer + aJac(nrgState,watState) = dVolHtCapBulk_dPsi0(iLayer) * mLayerTempPrime(jLayer) & + + (dt/mLayerDepth(jLayer))*(-dNrgFlux_dWatBelow(jLayer-1) + dNrgFlux_dWatAbove(jLayer)) + if(mLayerdTheta_dTk(jLayer) > verySmall)then ! ice is present + aJac(nrgState,watState) = -LH_fus*iden_water * dVolTot_dPsi0(iLayer) * cj & + - LH_fus*iden_water * mLayerMatricHeadPrime(iLayer) * d2VolTot_d2Psi0(iLayer) + aJac(nrgState,watState) ! dNrg/dMat (J m-3 m-1) -- dMat changes volumetric water, and hence ice content + endif + + ! - include derivatives of heat capacity w.r.t water fluxes for surrounding layers starting with layer above + if(iLayer>1)then + if(ixSoilOnlyNrg(iLayer-1)/=integerMissing) aJac(ixSoilOnlyNrg(iLayer-1),watState) = (dt/mLayerDepth(jLayer-1))*( dNrgFlux_dWatBelow(jLayer-1) ) + elseif(iLayer==1 .and. nSnowOnlyNrg>0)then !top soil layer and there is snow above + if(ixSnowOnlyNrg(nSnow)/=integerMissing) aJac(ixSnowOnlyNrg(nSnow),watState) = (dt/mLayerDepth(nSnow))*( dNrgFlux_dWatBelow(nSnow) ) + endif + + ! (cross-derivative terms for the layer below) + if(iLayer<nSoil)then + if(ixSoilOnlyHyd(iLayer+1)/=integerMissing) aJac(ixSoilOnlyNrg(iLayer+1),watState) = (dt/mLayerDepth(jLayer+1))*(-dNrgFlux_dWatAbove(jLayer ) ) + endif + + endif ! (if the water state for the current layer is within the state subset) + + ! - include terms for surface infiltration below surface + if(ixSoilOnlyHyd(1)/=integerMissing) aJac(ixSoilOnlyHyd(1),nrgState) = -(dt/mLayerDepth(1+nSnow))*dq_dNrgStateLayerSurfVec(iLayer) + aJac(ixSoilOnlyHyd(1),nrgState) + + end do ! (looping through soil layers) + + ! - include terms for surface infiltration above surface + if(nSnowOnlyHyd>0 .and. ixSnowOnlyNrg(nSnow)/=integerMissing)then + if(ixSoilOnlyHyd(1)/=integerMissing) aJac(ixSoilOnlyHyd(1),ixSnowOnlyNrg(nSnow)) = -(dt/mLayerDepth(1+nSnow))*dq_dNrgStateLayerSurfVec(0) + elseif(computeVegFlux .and. ixVegNrg/=integerMissing)then !ixTopHyd = ixSoilOnlyHyd(1) + if(ixTopHyd/=integerMissing) aJac(ixTopHyd,ixVegNrg) = -(dt/mLayerDepth(1+nSnow))*dq_dNrgStateLayerSurfVec(0) + aJac(ixTopHyd,ixVegNrg) + endif + + endif ! (if there are state variables for both water and energy in the soil domain) + + ! print the Jacobian + if(globalPrintFlag)then + print*, '** analytical Jacobian (full):' + write(*,'(a4,1x,100(i12,1x))') 'xCol', (iLayer, iLayer=min(iJac1,nState),min(iJac2,nState)) + do iLayer=min(iJac1,nState),min(iJac2,nState) + write(*,'(i4,1x,100(e12.5,1x))') iLayer, aJac(min(iJac1,nState):min(iJac2,nState),iLayer) + end do + endif + + !print*, '** analytical Jacobian (full):' + !write(*,'(a4,1x,100(i12,1x))') 'xCol', (iLayer, iLayer=1,size(aJac,2)) + !do iLayer=1,size(aJac,1) + ! write(*,'(i4,1x,100(e12.5,1x))') iLayer, aJac(iLayer,1:size(aJac,2)) + !end do + !print *, '--------------------------------------------------------------' + + ! *** + ! check + case default; err=20; message=trim(message)//'unable to identify option for the type of matrix'; return + + end select ! type of matrix + + if(any(isNan(aJac)))then + print *, '******************************* WE FOUND NAN IN JACOBIAN ************************************' + stop 1 + message=trim(message)//'we found NaN' + err=20; return + endif + + ! end association to variables in the data structures + end associate + + end subroutine computJacobSundials + + + ! ********************************************************************************************************** + ! public subroutine computJacobSetup: this sets up the inputs for the Jacobian computation for the IDA function + ! ********************************************************************************************************** + subroutine computJacobSetup(& + ! input: model control + cj, & ! intent(in): this scalar changes whenever the step size or method order changes + dt, & ! intent(in): time step + nSnow, & ! intent(in): number of snow layers + nSoil, & ! intent(in): number of soil layers + nLayers, & ! intent(in): total number of layers + ixMatrix, & ! intent(in): form of the Jacobian matrix + computeVegFlux, & ! intent(in): flag to indicate if we need to compute fluxes over vegetation + scalarSolution, & ! intent(in): flag to indicate the scalar solution + ! input: state vectors + stateVec, & ! intent(in): model state vector + stateVecPrime, & ! intent(in): derivative of model state vector + sMul, & ! intent(in): state vector multiplier (used in the residual calculations) + ! input: data structures + model_decisions, & ! intent(in): model decisions + mpar_data, & ! intent(in): model parameters + prog_data, & ! intent(in): model prognostic variables for a local HRU + ! input-output: data structures + indx_data, & ! intent(inout): index data + diag_data, & ! intent(inout): model diagnostic variables for a local HRU + deriv_data, & ! intent(inout): derivatives in model fluxes w.r.t. relevant state variables + ! input: baseflow + dBaseflow_dMatric, & ! intent(in): derivative in baseflow w.r.t. matric head (s-1) + ! output: flux and residual vectors + dMat, & ! intent(inout): diagonal of Jacobian Matrix + Jac, & ! intent(out): jacobian matrix + err,message) ! intent(out): error control + ! -------------------------------------------------------------------------------------------------------------------------------- + ! provide access to subroutines + USE getVectorzAddSundials_module, only:varExtractSundials + USE updateVarsSundials_module, only:updateVarsSundials ! update prognostic variables + implicit none + ! -------------------------------------------------------------------------------------------------------------------------------- + ! -------------------------------------------------------------------------------------------------------------------------------- + ! input: model control + real(rkind),intent(in) :: cj ! this scalar changes whenever the step size or method order changes + real(rkind),intent(in) :: dt ! time step + integer(i4b),intent(in) :: nSnow ! number of snow layers + integer(i4b),intent(in) :: nSoil ! number of soil layers + integer(i4b),intent(in) :: nLayers ! total number of layers + integer(i4b) :: ixMatrix ! form of matrix (band diagonal or full matrix) + logical(lgt),intent(in) :: computeVegFlux ! flag to indicate if computing fluxes over vegetation + logical(lgt),intent(in) :: scalarSolution ! flag to denote if implementing the scalar solution + ! input: state vectors + real(rkind),intent(in) :: stateVec(:) ! model state vector + real(rkind),intent(in) :: stateVecPrime(:) ! model state vector + real(qp),intent(in) :: sMul(:) ! NOTE: qp ! state vector multiplier (used in the residual calculations) + ! input: data structures + type(model_options),intent(in) :: model_decisions(:) ! model decisions + type(var_dlength), intent(in) :: mpar_data ! model parameters + type(var_dlength), intent(in) :: prog_data ! prognostic variables for a local HRU + ! output: data structures + type(var_ilength),intent(inout) :: indx_data ! indices defining model states and layers + type(var_dlength),intent(inout) :: diag_data ! diagnostic variables for a local HRU + type(var_dlength),intent(inout) :: deriv_data ! derivatives in model fluxes w.r.t. relevant state variables + ! input-output: baseflow + real(rkind),intent(in) :: dBaseflow_dMatric(:,:) ! derivative in baseflow w.r.t. matric head (s-1) + ! output: Jacobian + real(rkind), intent(inout) :: dMat(:) + real(rkind), intent(out) :: Jac(:,:) ! jacobian matrix + ! output: error control + integer(i4b),intent(out) :: err ! error code + character(*),intent(out) :: message ! error message + ! -------------------------------------------------------------------------------------------------------------------------------- + ! local variables + ! -------------------------------------------------------------------------------------------------------------------------------- + ! state variables + real(rkind) :: scalarCanairTempTrial ! trial value for temperature of the canopy air space (K) + real(rkind) :: scalarCanopyTempTrial ! trial value for temperature of the vegetation canopy (K) + real(rkind) :: scalarCanopyWatTrial ! trial value for liquid water storage in the canopy (kg m-2) + real(rkind),dimension(nLayers) :: mLayerTempTrial ! trial value for temperature of layers in the snow and soil domains (K) + real(rkind),dimension(nLayers) :: mLayerVolFracWatTrial ! trial value for volumetric fraction of total water (-) + real(rkind),dimension(nSoil) :: mLayerMatricHeadTrial ! trial value for total water matric potential (m) + real(rkind),dimension(nSoil) :: mLayerMatricHeadLiqTrial ! trial value for liquid water matric potential (m) + real(rkind) :: scalarAquiferStorageTrial ! trial value of storage of water in the aquifer (m) + ! diagnostic variables + real(rkind) :: scalarCanopyLiqTrial ! trial value for mass of liquid water on the vegetation canopy (kg m-2) + real(rkind) :: scalarCanopyIceTrial ! trial value for mass of ice on the vegetation canopy (kg m-2) + real(rkind),dimension(nLayers) :: mLayerVolFracLiqTrial ! trial value for volumetric fraction of liquid water (-) + real(rkind),dimension(nLayers) :: mLayerVolFracIceTrial ! trial value for volumetric fraction of ice (-) + ! derivative of state variables + real(rkind) :: scalarCanairTempPrime ! derivative value for temperature of the canopy air space (K) + real(rkind) :: scalarCanopyTempPrime ! derivative value for temperature of the vegetation canopy (K) + real(rkind) :: scalarCanopyWatPrime ! derivative value for liquid water storage in the canopy (kg m-2) + real(rkind),dimension(nLayers) :: mLayerTempPrime ! derivative value for temperature of layers in the snow and soil domains (K) + real(rkind),dimension(nLayers) :: mLayerVolFracWatPrime ! derivative value for volumetric fraction of total water (-) + real(rkind),dimension(nSoil) :: mLayerMatricHeadPrime ! derivative value for total water matric potential (m) + real(rkind),dimension(nSoil) :: mLayerMatricHeadLiqPrime ! derivative value for liquid water matric potential (m) + real(rkind) :: scalarAquiferStoragePrime ! derivative value of storage of water in the aquifer (m) + ! derivative of diagnostic variables + real(rkind) :: scalarCanopyLiqPrime ! derivative value for mass of liquid water on the vegetation canopy (kg m-2) + real(rkind) :: scalarCanopyIcePrime ! derivative value for mass of ice on the vegetation canopy (kg m-2) + real(rkind),dimension(nLayers) :: mLayerVolFracLiqPrime ! derivative value for volumetric fraction of liquid water (-) + real(rkind),dimension(nLayers) :: mLayerVolFracIcePrime ! derivative value for volumetric fraction of ice (-) + ! other local variables + character(LEN=256) :: cmessage ! error message of downwind routine + real(rkind) :: dt1 + + ! -------------------------------------------------------------------------------------------------------------------------------- + ! association to variables in the data structures + ! -------------------------------------------------------------------------------------------------------------------------------- + associate(& + ! model decisions + ixRichards => model_decisions(iLookDECISIONS%f_Richards)%iDecision ,& ! intent(in): [i4b] index of the form of Richards' equation + ! soil parameters + theta_sat => mpar_data%var(iLookPARAM%theta_sat)%dat ,& ! intent(in): [dp(:)] soil porosity (-) + specificStorage => mpar_data%var(iLookPARAM%specificStorage)%dat(1) ,& ! intent(in): [dp] specific storage coefficient (m-1) + ! decisions + ixGroundwater => model_decisions(iLookDECISIONS%groundwatr)%iDecision & + ) ! association to variables in the data structures + ! -------------------------------------------------------------------------------------------------------------------------------- + ! initialize error control + err=0; message="computJacobSetup/" + + ! these will need to be initialized as they do not have updated prognostic structures in Sundials + ! should all be set to previous values if splits, but for now operator splitting is not hooked up + scalarCanairTempPrime = realMissing + scalarCanopyTempPrime = realMissing + scalarCanopyWatPrime = realMissing + scalarCanopyLiqPrime = realMissing + scalarCanopyIcePrime = realMissing + mLayerTempPrime = realMissing + mLayerVolFracWatPrime = realMissing + mLayerVolFracLiqPrime = realMissing + mLayerVolFracIcePrime = realMissing + mLayerMatricHeadPrime = realMissing + mLayerMatricHeadLiqPrime = realMissing + scalarAquiferStoragePrime= realMissing + scalarCanairTempTrial = realMissing ! scalarCanairTemp prognostic not up to date + scalarCanopyTempTrial = realMissing ! scalarCanopyTempTrial prognostic not up to date + scalarCanopyWatTrial = realMissing ! scalarCanopyWat prognostic not up to date + scalarCanopyLiqTrial = realMissing ! scalarCanopyLiq prognostic not up to date + scalarCanopyIceTrial = realMissing ! scalarCanopyLIce prognostic not up to date + mLayerTempTrial = realMissing ! mLayerTemp prognostic not up to date + mLayerVolFracWatTrial = realMissing ! mLayerVolFracWat prognostic not up to date + mLayerVolFracLiqTrial = realMissing ! mLayerVolFracLiq prognostic not up to date + mLayerVolFracIceTrial = realMissing ! mLayerVolFracIce prognostic not up to date + mLayerMatricHeadTrial = realMissing ! mLayerMatricHead prognostic not up to date + mLayerMatricHeadLiqTrial = realMissing ! mLayerMatricHeadLiq prognostic not up to date + scalarAquiferStorageTrial= realMissing ! scalarAquiferStorage prognostic not up to date + + ! extract variables from the model state vector + call varExtractSundials(& + ! input + stateVec, & ! intent(in): model state vector (mixed units) + stateVecPrime, & ! intent(in): model state vector (mixed units) + diag_data, & ! intent(in): model diagnostic variables for a local HRU + prog_data, & ! intent(in): model prognostic variables for a local HRU + indx_data, & ! intent(in): indices defining model states and layers + ! output: variables for the vegetation canopy + scalarCanairTempTrial, & ! intent(inout): trial value of canopy air temperature (K) + scalarCanopyTempTrial, & ! intent(inout): trial value of canopy temperature (K) + scalarCanopyWatTrial, & ! intent(inout): trial value of canopy total water (kg m-2) + scalarCanopyLiqTrial, & ! intent(inout): trial value of canopy liquid water (kg m-2) + scalarCanairTempPrime, & ! intent(inout): derivative of canopy air temperature (K) + scalarCanopyTempPrime, & ! intent(inout): derivative of canopy temperature (K) + scalarCanopyWatPrime, & ! intent(inout): derivative of canopy total water (kg m-2) + scalarCanopyLiqPrime, & ! intent(inout): derivative of canopy liquid water (kg m-2) + ! output: variables for the snow-soil domain + mLayerTempTrial, & ! intent(inout): trial vector of layer temperature (K) + mLayerVolFracWatTrial, & ! intent(inout): trial vector of volumetric total water content (-) + mLayerVolFracLiqTrial, & ! intent(inout): trial vector of volumetric liquid water content (-) + mLayerMatricHeadTrial, & ! intent(inout): trial vector of total water matric potential (m) + mLayerMatricHeadLiqTrial, & ! intent(inout): trial vector of liquid water matric potential (m) + mLayerTempPrime, & ! intent(inout): derivative of layer temperature (K) + mLayerVolFracWatPrime, & ! intent(inout): derivative of volumetric total water content (-) + mLayerVolFracLiqPrime, & ! intent(inout): derivative of volumetric liquid water content (-) + mLayerMatricHeadPrime, & ! intent(inout): derivative of total water matric potential (m) + mLayerMatricHeadLiqPrime, & ! intent(inout): derivative of liquid water matric potential (m) + ! output: variables for the aquifer + scalarAquiferStorageTrial,& ! intent(inout): trial value of storage of water in the aquifer (m) + scalarAquiferStoragePrime,& ! intent(inout): derivative of storage of water in the aquifer (m) + ! output: error control + err,cmessage) ! intent(out): error control + if(err/=0)then; message=trim(message)//trim(cmessage); return; end if ! (check for errors) + + call updateVarsSundials(& + ! input + dt, & ! intent(in): time step + .true., & ! intent(in): logical flag if computing Jacobian for Sundials solver + .false., & ! intent(in): logical flag to adjust temperature to account for the energy used in melt+freeze + mpar_data, & ! intent(in): model parameters for a local HRU + indx_data, & ! intent(in): indices defining model states and layers + prog_data, & ! intent(in): model prognostic variables for a local HRU + mLayerVolFracWatTrial, & ! intent(in): use current vector for prev vector of volumetric total water content (-) + mLayerMatricHeadTrial, & ! intent(in): use current vector for prev vector of total water matric potential (m) + diag_data, & ! intent(inout): model diagnostic variables for a local HRU + deriv_data, & ! intent(inout): derivatives in model fluxes w.r.t. relevant state variables + ! output: variables for the vegetation canopy + scalarCanopyTempTrial, & ! intent(inout): trial value of canopy temperature (K) + scalarCanopyWatTrial, & ! intent(inout): trial value of canopy total water (kg m-2) + scalarCanopyLiqTrial, & ! intent(inout): trial value of canopy liquid water (kg m-2) + scalarCanopyIceTrial, & ! intent(inout): trial value of canopy ice content (kg m-2) + scalarCanopyTempPrime, & ! intent(inout): trial value of canopy temperature (K) + scalarCanopyWatPrime, & ! intent(inout): trial value of canopy total water (kg m-2) + scalarCanopyLiqPrime, & ! intent(inout): trial value of canopy liquid water (kg m-2) + scalarCanopyIcePrime, & ! intent(inout): trial value of canopy ice content (kg m-2 + ! output: variables for the snow-soil domain + mLayerTempTrial, & ! intent(inout): trial vector of layer temperature (K) + mLayerVolFracWatTrial, & ! intent(inout): trial vector of volumetric total water content (-) + mLayerVolFracLiqTrial, & ! intent(inout): trial vector of volumetric liquid water content (-) + mLayerVolFracIceTrial, & ! intent(inout): trial vector of volumetric ice water content (-) + mLayerMatricHeadTrial, & ! intent(inout): trial vector of total water matric potential (m) + mLayerMatricHeadLiqTrial, & ! intent(inout): trial vector of liquid water matric potential (m) + mLayerTempPrime, & ! + mLayerVolFracWatPrime, & ! intent(inout): Prime vector of volumetric total water content (-) + mLayerVolFracLiqPrime, & ! intent(inout): Prime vector of volumetric liquid water content (-) + mLayerVolFracIcePrime, & ! + mLayerMatricHeadPrime, & ! intent(inout): Prime vector of total water matric potential (m) + mLayerMatricHeadLiqPrime, & ! intent(inout): Prime vector of liquid water matric potential (m) + ! output: error control + err,cmessage) ! intent(out): error control + if(err/=0)then; message=trim(message)//trim(cmessage); return; end if ! (check for errors) + + ! ----- + ! * compute the Jacobian matrix... + ! -------------------------------- + + ! compute the analytical Jacobian matrix + ! NOTE: The derivatives were computed in the previous call to computFlux + ! This occurred either at the call to eval8summaSundials at the start of systemSolvSundials + ! or in the call to eval8summaSundials in the previous iteration + dt1 = 1._qp + call computJacobSundials(& + ! input: model control + cj, & ! intent(in): this scalar changes whenever the step size or method order changes + dt1, & ! intent(in): length of the time step (seconds) + nSnow, & ! intent(in): number of snow layers + nSoil, & ! intent(in): number of soil layers + nLayers, & ! intent(in): total number of layers + computeVegFlux, & ! intent(in): flag to indicate if we need to compute fluxes over vegetation + (ixGroundwater==qbaseTopmodel), & ! intent(in): flag to indicate if we need to compute baseflow + ixMatrix, & ! intent(in): form of the Jacobian matrix + specificStorage, & ! intent(in): specific storage coefficient (m-1) + theta_sat, & ! intent(in): soil porosity (-) + ixRichards, & ! intent(in): choice of option for Richards' equation + ! input: data structures + indx_data, & ! intent(in): index data + prog_data, & ! intent(in): model prognostic variables for a local HRU + diag_data, & ! intent(in): model diagnostic variables for a local HRU + deriv_data, & ! intent(in): derivatives in model fluxes w.r.t. relevant state variables + dBaseflow_dMatric, & ! intent(in): derivative in baseflow w.r.t. matric head (s-1) + ! input: state variables + mLayerTempPrime, & ! intent(in) + mLayerMatricHeadPrime, & ! intent(in) + mLayerMatricHeadLiqPrime, & ! intent(in) + mLayerVolFracWatPrime, & ! intent(in) + scalarCanopyTempPrime, & ! intent(in) derivative value for temperature of the vegetation canopy (K) + scalarCanopyWatPrime, & ! intent(in) + ! input-output: Jacobian and its diagonal + dMat, & ! intent(inout): diagonal of the Jacobian matrix + Jac, & ! intent(out): Jacobian matrix + ! output: error control + err,cmessage) ! intent(out): error code and error message + if(err/=0)then; message=trim(message)//trim(cmessage); return; end if ! (check for errors) + + ! end association with the information in the data structures + end associate + + end subroutine computJacobSetup + + + ! ********************************************************************************************************** + ! public function computJacob4IDA: the interface to compute the Jacobian matrix dF/dy + c dF/dy' for IDA solver + ! ********************************************************************************************************** + ! Return values: + ! 0 = success, + ! 1 = recoverable error, + ! -1 = non-recoverable error + ! ---------------------------------------------------------------- + integer(c_int) function computJacob4IDA(t, cj, sunvec_y, sunvec_yp, sunvec_r, & + sunmat_J, user_data, sunvec_temp1, sunvec_temp2, sunvec_temp3) & + result(ierr) bind(C,name='computJacob4IDA') + + !======= Inclusions =========== + use, intrinsic :: iso_c_binding + use fsundials_nvector_mod + use fsundials_matrix_mod + use fnvector_serial_mod + use fsunmatrix_dense_mod + use nrtype + use type4IDA + !======= Declarations ========= + implicit none + + ! calling variables + real(rkind), value :: t ! current time + real(rkind), value :: cj ! step size scaling factor + type(N_Vector) :: sunvec_y ! solution N_Vector + type(N_Vector) :: sunvec_yp ! derivative N_Vector + type(N_Vector) :: sunvec_r ! residual N_Vector + type(SUNMatrix) :: sunmat_J ! Jacobian SUNMatrix + type(c_ptr), value :: user_data ! user-defined data + type(N_Vector) :: sunvec_temp1 ! temporary N_Vector + type(N_Vector) :: sunvec_temp2 ! temporary N_Vector + type(N_Vector) :: sunvec_temp3 ! temporary N_Vector + + ! pointers to data in SUNDIALS vectors + real(rkind), pointer :: stateVec(:) ! state vector + real(rkind), pointer :: stateVecPrime(:)! derivative of the state vector + real(rkind), pointer :: rVec(:) ! residual vector + real(rkind), pointer :: Jac(:,:) ! Jacobian matrix + type(eqnsData), pointer :: eqns_data ! equations data + + + + !======= Internals ============ + + ! get equations data from user-defined data + call c_f_pointer(user_data, eqns_data) + + + ! get data arrays from SUNDIALS vectors + stateVec(1:eqns_data%nState) => FN_VGetArrayPointer(sunvec_y) + stateVecPrime(1:eqns_data%nState) => FN_VGetArrayPointer(sunvec_yp) + rVec(1:eqns_data%nState) => FN_VGetArrayPointer(sunvec_r) + Jac(1:eqns_data%nState, 1:eqns_data%nState) => FSUNDenseMatrix_Data(sunmat_J) + + ! compute Jacobian matrix + call computJacobSetup(& + ! input: model control + cj, & ! intent(in): this scalar changes whenever the step size or method order changes + eqns_data%dt, & ! intent(in): data step + eqns_data%nSnow, & ! intent(in): number of snow layers + eqns_data%nSoil, & ! intent(in): number of soil layers + eqns_data%nLayers, & ! intent(in): number of layers + eqns_data%ixMatrix, & ! intent(in): type of matrix (dense or banded) + eqns_data%computeVegFlux, & ! intent(in): flag to indicate if we need to compute fluxes over vegetation + eqns_data%scalarSolution, & ! intent(in): flag to indicate the scalar solution + ! input: state vectors + stateVec, & ! intent(in): model state vector + stateVecPrime, & ! intent(in): model state vector + eqns_data%sMul, & ! intent(in): state vector multiplier (used in the residual calculations) + ! input: data structures + model_decisions, & ! intent(in): model decisions + eqns_data%mpar_data, & ! intent(in): model parameters + eqns_data%prog_data, & ! intent(in): model prognostic variables for a local HRU + ! input-output: data structures + eqns_data%indx_data, & ! intent(inou): index data + eqns_data%diag_data, & ! intent(inout): model diagnostic variables for a local HRU + eqns_data%deriv_data, & ! intent(inout): derivatives in model fluxes w.r.t. relevant state variables + ! input: baseflow + eqns_data%dBaseflow_dMatric, & ! intent(in): derivative in baseflow w.r.t. matric head (s-1) + ! output + eqns_data%dMat, & ! intetn(inout): diagonal of the Jacobian matrix + Jac, & ! intent(out): Jacobian matrix + eqns_data%err,eqns_data%message) ! intent(out): error control + + if(eqns_data%err > 0)then; eqns_data%message=trim(eqns_data%message); ierr=-1; return; endif + if(eqns_data%err < 0)then; eqns_data%message=trim(eqns_data%message); ierr=1; return; endif + + ! return success + ierr = 0 + return + + + + end function computJacob4IDA + + + ! ********************************************************************************************************** + ! private function: get the off-diagonal index in the band-diagonal matrix + ! ********************************************************************************************************** + function ixOffDiag(jState,iState) + implicit none + integer(i4b),intent(in) :: jState ! off-diagonal state + integer(i4b),intent(in) :: iState ! diagonal state + integer(i4b) :: ixOffDiag ! off-diagonal index in gthe band-diagonal matrix + ixOffDiag = ixBandOffset + jState - iState + end function ixOffDiag + + end module computJacobSundials_module + \ No newline at end of file diff --git a/build/source/engine/sundials/computResidDAE.f90 b/build/source/engine/sundials/computResidDAE.f90 deleted file mode 100644 index a87b2e7..0000000 --- a/build/source/engine/sundials/computResidDAE.f90 +++ /dev/null @@ -1,346 +0,0 @@ - - -module computResidDAE_module - -! data types -USE nrtype - -! derived types to define the data structures -USE data_types,only:& - var_ilength, & ! data vector with variable length dimension (i4b) - var_dlength ! data vector with variable length dimension (rkind) - -! named variables -USE var_lookup,only:iLookPROG ! named variables for structure elements -USE var_lookup,only:iLookDIAG ! named variables for structure elements -USE var_lookup,only:iLookFLUX ! named variables for structure elements -USE var_lookup,only:iLookINDEX ! named variables for structure elements - -! access the global print flag -USE globalData,only:globalPrintFlag - -! access missing values -USE globalData,only:integerMissing ! missing integer -USE globalData,only:realMissing ! missing real number - -! define access to state variables to print -USE globalData,only: iJac1 ! first layer of the Jacobian to print -USE globalData,only: iJac2 ! last layer of the Jacobian to print - -! domain types -USE globalData,only:iname_veg ! named variables for vegetation -USE globalData,only:iname_snow ! named variables for snow -USE globalData,only:iname_soil ! named variables for soil - -! named variables to describe the state variable type -USE globalData,only:iname_nrgCanair ! named variable defining the energy of the canopy air space -USE globalData,only:iname_nrgCanopy ! named variable defining the energy of the vegetation canopy -USE globalData,only:iname_watCanopy ! named variable defining the mass of water on the vegetation canopy -USE globalData,only:iname_nrgLayer ! named variable defining the energy state variable for snow+soil layers -USE globalData,only:iname_watLayer ! named variable defining the total water state variable for snow+soil layers -USE globalData,only:iname_liqLayer ! named variable defining the liquid water state variable for snow+soil layers -USE globalData,only:iname_matLayer ! named variable defining the matric head state variable for soil layers -USE globalData,only:iname_lmpLayer ! named variable defining the liquid matric potential state variable for soil layers - -! constants -USE multiconst,only:& - LH_fus, & ! latent heat of fusion (J kg-1) - iden_ice, & ! intrinsic density of ice (kg m-3) - iden_water ! intrinsic density of liquid water (kg m-3) -! privacy -implicit none -private::printResidDAE -public::computResidDAE -contains - - ! ********************************************************************************************************** - ! public subroutine computResidDAE: compute the residual vector - ! ********************************************************************************************************** - subroutine computResidDAE(& - ! input: model control - nSnow, & ! intent(in): number of snow layers - nSoil, & ! intent(in): number of soil layers - nLayers, & ! intent(in): total number of layers - ! input: flux vectors - sMul, & ! intent(in): state vector multiplier (used in the residual calculations) - fVec, & ! intent(in): flux vector - ! input: state variables (already disaggregated into scalars and vectors) - scalarCanopyTempTrial, & ! intent(in): - mLayerTempTrial, & ! intent(in) - scalarCanairTempPrime, & ! intent(in): trial value for the temperature of the canopy air space (K) - scalarCanopyTempPrime, & ! intent(in): trial value for the temperature of the vegetation canopy (K) - scalarCanopyWatPrime, & ! - mLayerTempPrime, & ! intent(in): trial value for the temperature of each snow and soil layer (K) water content - scalarAquiferStoragePrime, & ! intent(in): trial value of storage of water in the aquifer (m) - ! input: diagnostic variables defining the liquid water and ice content (function of state variables) - scalarCanopyIcePrime, & ! intent(in): trial value for the ice on the vegetation canopy (kg m-2) - scalarCanopyLiqPrime, & ! intent(in): - mLayerVolFracIcePrime, & ! intent(in): trial value for the volumetric ice in each snow and soil layer (-) - mLayerVolFracWatPrime, & - mLayerVolFracLiqPrime, & - scalarCanopyCmTrial, & - mLayerCmTrial, & ! intent(in) - ! input: data structures - prog_data, & ! intent(in): model prognostic variables for a local HRU - diag_data, & ! intent(in): model diagnostic variables for a local HRU - flux_data, & ! intent(in): model fluxes for a local HRU - indx_data, & ! intent(in): index data - ! output - rAdd, & ! intent(out): additional (sink) terms on the RHS of the state equation - rVec, & ! intent(out): residual vector - err,message) ! intent(out): error control - ! -------------------------------------------------------------------------------------------------------------------------------- - implicit none - ! input: model control - integer(i4b),intent(in) :: nSnow ! number of snow layers - integer(i4b),intent(in) :: nSoil ! number of soil layers - integer(i4b),intent(in) :: nLayers ! total number of layers in the snow+soil domain - ! input: flux vectors - real(qp),intent(in) :: sMul(:) ! NOTE: qp ! state vector multiplier (used in the residual calculations) - real(rkind),intent(in) :: fVec(:) ! flux vector - ! input: state variables (already disaggregated into scalars and vectors) - real(rkind),intent(in) :: scalarCanopyTempTrial - real(rkind),intent(in) :: mLayerTempTrial(:) - real(rkind),intent(in) :: scalarCanairTempPrime ! trial value for temperature of the canopy air space (K) - real(rkind),intent(in) :: scalarCanopyTempPrime ! trial value for temperature of the vegetation canopy (K) - real(rkind),intent(in) :: scalarCanopyWatPrime ! derivative value for liquid water storage in the canopy (kg m-2) - real(rkind),intent(in) :: mLayerTempPrime(:) ! trial value for temperature of each snow/soil layer (K) content - real(rkind),intent(in) :: scalarAquiferStoragePrime ! trial value of aquifer storage (m) - ! input: diagnostic variables defining the liquid water and ice content (function of state variables) - real(rkind),intent(in) :: scalarCanopyIcePrime ! trial value for mass of ice on the vegetation canopy (kg m-2) - real(rkind),intent(in) :: scalarCanopyLiqPrime - real(rkind),intent(in) :: mLayerVolFracIcePrime(:) ! trial value for volumetric fraction of ice (-) - real(rkind),intent(in) :: mLayerVolFracLiqPrime(:) - real(rkind),intent(in) :: mLayerVolFracWatPrime(:) - real(qp),intent(in) :: scalarCanopyCmTrial - real(qp),intent(in) :: mLayerCmTrial(:) - ! input: data structures - type(var_dlength),intent(in) :: prog_data ! prognostic variables for a local HRU - type(var_dlength),intent(in) :: diag_data ! diagnostic variables for a local HRU - type(var_dlength),intent(in) :: flux_data ! model fluxes for a local HRU - type(var_ilength),intent(in) :: indx_data ! indices defining model states and layers - ! output - real(rkind),intent(out) :: rAdd(:) ! additional (sink) terms on the RHS of the state equation - real(qp),intent(out) :: rVec(:) ! NOTE: qp ! residual vector - integer(i4b),intent(out) :: err ! error code - character(*),intent(out) :: message ! error message - ! -------------------------------------------------------------------------------------------------------------------------------- - ! local variables - ! -------------------------------------------------------------------------------------------------------------------------------- - integer(i4b) :: iLayer ! index of layer within the snow+soil domain - integer(i4b),parameter :: ixVegVolume=1 ! index of the desired vegetation control volumne (currently only one veg layer) - real(rkind),dimension(nLayers) :: mLayerVolFracHydPrime ! vector of volumetric water content (-), either liquid water content or total water content - real(rkind) :: scalarCanopyHydPrime ! trial value for canopy water (kg m-2), either liquid water content or total water content - ! -------------------------------------------------------------------------------------------------------------------------------- - ! -------------------------------------------------------------------------------------------------------------------------------- - ! link to the necessary variables for the residual computations - associate(& - ! canopy and layer depth - canopyDepth => diag_data%var(iLookDIAG%scalarCanopyDepth)%dat(1) ,& ! intent(in): [dp] canopy depth (m) - mLayerDepth => prog_data%var(iLookPROG%mLayerDepth)%dat ,& ! intent(in): [dp(:)] depth of each layer in the snow-soil sub-domain (m) - ! model fluxes (sink terms in the soil domain) - mLayerTranspire => flux_data%var(iLookFLUX%mLayerTranspire)%dat ,& ! intent(in): [dp] transpiration loss from each soil layer (m s-1) - mLayerBaseflow => flux_data%var(iLookFLUX%mLayerBaseflow)%dat ,& ! intent(in): [dp(:)] baseflow from each soil layer (m s-1) - mLayerCompress => diag_data%var(iLookDIAG%mLayerCompress)%dat ,& ! intent(in): [dp(:)] change in storage associated with compression of the soil matrix (-) - ! number of state variables of a specific type - nSnowSoilNrg => indx_data%var(iLookINDEX%nSnowSoilNrg )%dat(1) ,& ! intent(in): [i4b] number of energy state variables in the snow+soil domain - nSnowSoilHyd => indx_data%var(iLookINDEX%nSnowSoilHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the snow+soil domain - nSoilOnlyHyd => indx_data%var(iLookINDEX%nSoilOnlyHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the soil domain - ! model indices - ixCasNrg => indx_data%var(iLookINDEX%ixCasNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy air space energy state variable - ixVegNrg => indx_data%var(iLookINDEX%ixVegNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy energy state variable - ixVegHyd => indx_data%var(iLookINDEX%ixVegHyd)%dat(1) ,& ! intent(in): [i4b] index of canopy hydrology state variable (mass) - ixAqWat => indx_data%var(iLookINDEX%ixAqWat)%dat(1) ,& ! intent(in): [i4b] index of water storage in the aquifer - ixSnowSoilNrg => indx_data%var(iLookINDEX%ixSnowSoilNrg)%dat ,& ! intent(in): [i4b(:)] indices for energy states in the snow+soil subdomain - ixSnowSoilHyd => indx_data%var(iLookINDEX%ixSnowSoilHyd)%dat ,& ! intent(in): [i4b(:)] indices for hydrology states in the snow+soil subdomain - ixSoilOnlyHyd => indx_data%var(iLookINDEX%ixSoilOnlyHyd)%dat ,& ! intent(in): [i4b(:)] indices for hydrology states in the soil subdomain - ixStateType => indx_data%var(iLookINDEX%ixStateType)%dat ,& ! intent(in): [i4b(:)] indices defining the type of the state (iname_nrgLayer...) - ixHydCanopy => indx_data%var(iLookINDEX%ixHydCanopy)%dat ,& ! intent(in): [i4b(:)] index of the hydrology states in the canopy domain - ixHydType => indx_data%var(iLookINDEX%ixHydType)%dat ,& ! intent(in): [i4b(:)] named variables defining the type of hydrology states in snow+soil domain - layerType => indx_data%var(iLookINDEX%layerType)%dat & ! intent(in): [i4b(:)] named variables defining the type of layer in snow+soil domain - ) ! association to necessary variables for the residual computations - ! -------------------------------------------------------------------------------------------------------------------------------- - ! initialize error control - err=0; message="computResidDAE/" - - ! --- - ! * compute sink terms... - ! ----------------------- - - ! intialize additional terms on the RHS as zero - rAdd(:) = 0._rkind - - ! compute energy associated with melt freeze for the vegetation canopy - if(ixVegNrg/=integerMissing) rAdd(ixVegNrg) = rAdd(ixVegNrg) + LH_fus*scalarCanopyIcePrime/canopyDepth ! energy associated with melt/freeze (J m-3) - ! compute energy associated with melt/freeze for snow - ! NOTE: allow expansion of ice during melt-freeze for snow; deny expansion of ice during melt-freeze for soil - if(nSnowSoilNrg>0)then - do concurrent (iLayer=1:nLayers,ixSnowSoilNrg(iLayer)/=integerMissing) ! (loop through non-missing energy state variables in the snow+soil domain) - select case( layerType(iLayer) ) - case(iname_snow); rAdd( ixSnowSoilNrg(iLayer) ) = rAdd( ixSnowSoilNrg(iLayer) ) + LH_fus*iden_ice * mLayerVolFracIcePrime(iLayer) - case(iname_soil); rAdd( ixSnowSoilNrg(iLayer) ) = rAdd( ixSnowSoilNrg(iLayer) ) + LH_fus*iden_water * mLayerVolFracIcePrime(iLayer) - end select - end do ! looping through non-missing energy state variables in the snow+soil domain - endif - - ! sink terms soil hydrology (-) - ! NOTE 1: state variable is volumetric water content, so melt-freeze is not included - ! NOTE 2: ground evaporation was already included in the flux at the upper boundary - ! NOTE 3: rAdd(ixSnowOnlyWat)=0, and is defined in the initialization above - ! NOTE 4: same sink terms for matric head and liquid matric potential - if(nSoilOnlyHyd>0)then - do concurrent (iLayer=1:nSoil,ixSoilOnlyHyd(iLayer)/=integerMissing) ! (loop through non-missing hydrology state variables in the snow+soil domain) - rAdd( ixSoilOnlyHyd(iLayer) ) = rAdd( ixSoilOnlyHyd(iLayer) ) + (mLayerTranspire(iLayer) - mLayerBaseflow(iLayer) )/mLayerDepth(iLayer+nSnow) - mLayerCompress(iLayer) - end do ! looping through non-missing energy state variables in the snow+soil domain - endif - - ! --- - ! * compute the residual vector... - ! -------------------------------- - - ! compute the residual vector for the vegetation canopy - ! NOTE: sMul(ixVegHyd) = 1, but include as it converts all variables to quadruple precision - ! --> energy balance - if(ixCasNrg/=integerMissing) rVec(ixCasNrg) = sMul(ixCasNrg)*scalarCanairTempPrime - ( fVec(ixCasNrg) + rAdd(ixCasNrg) ) - if(ixVegNrg/=integerMissing) rVec(ixVegNrg) = sMul(ixVegNrg) * scalarCanopyTempPrime + scalarCanopyCmTrial * scalarCanopyWatPrime/canopyDepth - ( fVec(ixVegNrg) + rAdd(ixVegNrg) ) - ! --> mass balance - if(ixVegHyd/=integerMissing)then - scalarCanopyHydPrime = merge(scalarCanopyWatPrime, scalarCanopyLiqPrime, (ixStateType( ixHydCanopy(ixVegVolume) )==iname_watCanopy) ) - rVec(ixVegHyd) = sMul(ixVegHyd)*scalarCanopyHydPrime - ( fVec(ixVegHyd) + rAdd(ixVegHyd) ) - endif - - ! compute the residual vector for the snow and soil sub-domains for energy - if(nSnowSoilNrg>0)then - do concurrent (iLayer=1:nLayers,ixSnowSoilNrg(iLayer)/=integerMissing) ! (loop through non-missing energy state variables in the snow+soil domain) - rVec( ixSnowSoilNrg(iLayer) ) = sMul( ixSnowSoilNrg(iLayer) ) * mLayerTempPrime(iLayer) + mLayerCmTrial(iLayer) * mLayerVolFracWatPrime(iLayer) - ( fVec( ixSnowSoilNrg(iLayer) ) + rAdd( ixSnowSoilNrg(iLayer) ) ) - end do ! looping through non-missing energy state variables in the snow+soil domain - endif - - ! compute the residual vector for the snow and soil sub-domains for hydrology - ! NOTE: residual depends on choice of state variable - if(nSnowSoilHyd>0)then - do concurrent (iLayer=1:nLayers,ixSnowSoilHyd(iLayer)/=integerMissing) ! (loop through non-missing hydrology state variables in the snow+soil domain) - ! (get the correct state variable) - mLayerVolFracHydPrime(iLayer) = merge(mLayerVolFracWatPrime(iLayer), mLayerVolFracLiqPrime(iLayer), (ixHydType(iLayer)==iname_watLayer .or. ixHydType(iLayer)==iname_matLayer) ) - ! (compute the residual) - rVec( ixSnowSoilHyd(iLayer) ) = mLayerVolFracHydPrime(iLayer) - ( fVec( ixSnowSoilHyd(iLayer) ) + rAdd( ixSnowSoilHyd(iLayer) ) ) - end do ! looping through non-missing energy state variables in the snow+soil domain - endif - - ! compute the residual vector for the aquifer - if(ixAqWat/=integerMissing) rVec(ixAqWat) = sMul(ixAqWat)*scalarAquiferStoragePrime - ( fVec(ixAqWat) + rAdd(ixAqWat) ) - - ! print result - if(globalPrintFlag)then - write(*,'(a,1x,100(e12.5,1x))') 'rVec = ', rVec(min(iJac1,size(rVec)):min(iJac2,size(rVec))) - write(*,'(a,1x,100(e12.5,1x))') 'fVec = ', fVec(min(iJac1,size(rVec)):min(iJac2,size(rVec))) - !print*, 'PAUSE:'; read(*,*) - endif - - !call printResidDAE(nSnow,nSoil,nLayers,indx_data,rAdd,rVec) - - ! check - if(any(isNan(rVec)))then - call printResidDAE(nSnow,nSoil,nLayers,indx_data,rAdd,rVec) - message=trim(message)//'we found NaN' - err=20; return - endif - - ! end association with the necessary variabiles for the residual calculations - end associate - - end subroutine computResidDAE - - ! ********************************************************************************************************** - ! private subroutine printResidDAE: print the residual vector mainly for debugging - ! ********************************************************************************************************** - subroutine printResidDAE( & - ! input: model control - nSnow, & ! intent(in): number of snow layers - nSoil, & ! intent(in): number of soil layers - nLayers, & ! intent(in): total number of layers - ! input: data structures - indx_data, & ! intent(in): index data - ! output - rAdd, & ! intent(out): additional (sink) terms on the RHS of the state equation - rVec) ! intent(out): residual vector - -! -------------------------------------------------------------------------------------------------------------------------------- -implicit none -! input: model control -integer(i4b),intent(in) :: nSnow ! number of snow layers -integer(i4b),intent(in) :: nSoil ! number of soil layers -integer(i4b),intent(in) :: nLayers ! total number of layers in the snow+soil domain -type(var_ilength),intent(in) :: indx_data ! indices defining model states and layers -! output -real(rkind),intent(in) :: rAdd(:) ! additional (sink) terms on the RHS of the state equation -real(qp),intent(in) :: rVec(:) ! NOTE: qp ! residual vector -! -------------------------------------------------------------------------------------------------------------------------------- -! local variables -! -------------------------------------------------------------------------------------------------------------------------------- -integer(i4b) :: iLayer ! index of layer within the snow+soil domain -! -------------------------------------------------------------------------------------------------------------------------------- -! link to the necessary variables for the residual computations -associate(& -! number of state variables of a specific type -nSnowSoilNrg => indx_data%var(iLookINDEX%nSnowSoilNrg )%dat(1) ,& ! intent(in): [i4b] number of energy state variables in the snow+soil domain -nSnowSoilHyd => indx_data%var(iLookINDEX%nSnowSoilHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the snow+soil domain -nSoilOnlyHyd => indx_data%var(iLookINDEX%nSoilOnlyHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the soil domain -! model indices -ixCasNrg => indx_data%var(iLookINDEX%ixCasNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy air space energy state variable -ixVegNrg => indx_data%var(iLookINDEX%ixVegNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy energy state variable -ixVegHyd => indx_data%var(iLookINDEX%ixVegHyd)%dat(1) ,& ! intent(in): [i4b] index of canopy hydrology state variable (mass) -ixAqWat => indx_data%var(iLookINDEX%ixAqWat)%dat(1) ,& ! intent(in): [i4b] index of water storage in the aquifer -ixSnowSoilNrg => indx_data%var(iLookINDEX%ixSnowSoilNrg)%dat ,& ! intent(in): [i4b(:)] indices for energy states in the snow+soil subdomain -ixSnowSoilHyd => indx_data%var(iLookINDEX%ixSnowSoilHyd)%dat ,& ! intent(in): [i4b(:)] indices for hydrology states in the snow+soil subdomain -ixSoilOnlyHyd => indx_data%var(iLookINDEX%ixSoilOnlyHyd)%dat ,& ! intent(in): [i4b(:)] indices for hydrology states in the soil subdomain -ixStateType => indx_data%var(iLookINDEX%ixStateType)%dat ,& ! intent(in): [i4b(:)] indices defining the type of the state (iname_nrgLayer...) -ixHydCanopy => indx_data%var(iLookINDEX%ixHydCanopy)%dat ,& ! intent(in): [i4b(:)] index of the hydrology states in the canopy domain -ixHydType => indx_data%var(iLookINDEX%ixHydType)%dat ,& ! intent(in): [i4b(:)] named variables defining the type of hydrology states in snow+soil domain -layerType => indx_data%var(iLookINDEX%layerType)%dat & ! intent(in): [i4b(:)] named variables defining the type of layer in snow+soil domain -) ! association to necessary variables for the residual computations -! -------------------------------------------------------------------------------------------------------------------------------- - -if(ixVegNrg/=integerMissing) print *, 'rAdd(ixVegNrg) = ', rAdd(ixVegNrg) - -if(nSnowSoilNrg>0)then - do concurrent (iLayer=1:nLayers,ixSnowSoilNrg(iLayer)/=integerMissing) - select case( layerType(iLayer) ) - case(iname_snow) - print *, 'rAdd( ixSnowSoilNrg(iLayer) ) = ', rAdd( ixSnowSoilNrg(iLayer) ) - case(iname_soil); print *, 'rAdd( ixSnowSoilNrg(iLayer) ) = ', rAdd( ixSnowSoilNrg(iLayer) ) - end select - end do -endif - -if(nSoilOnlyHyd>0)then - do concurrent (iLayer=1:nSoil,ixSoilOnlyHyd(iLayer)/=integerMissing) - print *, 'rAdd( ixSoilOnlyHyd(iLayer) ) = ', rAdd( ixSoilOnlyHyd(iLayer) ) - end do -endif - -if(ixCasNrg/=integerMissing) print *, 'rVec(ixCasNrg) = ', rVec(ixCasNrg) -if(ixVegNrg/=integerMissing) print *, 'rVec(ixVegNrg) = ', rVec(ixVegNrg) -if(ixVegHyd/=integerMissing)then - print *, 'rVec(ixVegHyd) = ', rVec(ixVegHyd) -endif - -if(nSnowSoilNrg>0)then - do concurrent (iLayer=1:nLayers,ixSnowSoilNrg(iLayer)/=integerMissing) - print *, 'rVec( ixSnowSoilNrg(iLayer) ) = ', rVec( ixSnowSoilNrg(iLayer) ) - end do -endif - -if(nSnowSoilHyd>0)then - do concurrent (iLayer=1:nLayers,ixSnowSoilHyd(iLayer)/=integerMissing) - print *, 'rVec( ixSnowSoilHyd(iLayer) ) = ', rVec( ixSnowSoilHyd(iLayer) ) - end do -endif - -if(ixAqWat/=integerMissing) print *, ' rVec(ixAqWat) = ', rVec(ixAqWat) - -end associate - -end subroutine printResidDAE - -end module computResidDAE_module diff --git a/build/source/engine/sundials/computResidSundials.f90 b/build/source/engine/sundials/computResidSundials.f90 new file mode 100644 index 0000000..503f082 --- /dev/null +++ b/build/source/engine/sundials/computResidSundials.f90 @@ -0,0 +1,347 @@ + + +module computResidSundials_module + + ! data types + USE nrtype + + ! derived types to define the data structures + USE data_types,only:& + var_ilength, & ! data vector with variable length dimension (i4b) + var_dlength ! data vector with variable length dimension (rkind) + + ! named variables + USE var_lookup,only:iLookPROG ! named variables for structure elements + USE var_lookup,only:iLookDIAG ! named variables for structure elements + USE var_lookup,only:iLookFLUX ! named variables for structure elements + USE var_lookup,only:iLookINDEX ! named variables for structure elements + + ! access the global print flag + USE globalData,only:globalPrintFlag + + ! access missing values + USE globalData,only:integerMissing ! missing integer + USE globalData,only:realMissing ! missing real number + + ! define access to state variables to print + USE globalData,only: iJac1 ! first layer of the Jacobian to print + USE globalData,only: iJac2 ! last layer of the Jacobian to print + + ! domain types + USE globalData,only:iname_veg ! named variables for vegetation + USE globalData,only:iname_snow ! named variables for snow + USE globalData,only:iname_soil ! named variables for soil + + ! named variables to describe the state variable type + USE globalData,only:iname_nrgCanair ! named variable defining the energy of the canopy air space + USE globalData,only:iname_nrgCanopy ! named variable defining the energy of the vegetation canopy + USE globalData,only:iname_watCanopy ! named variable defining the mass of water on the vegetation canopy + USE globalData,only:iname_nrgLayer ! named variable defining the energy state variable for snow+soil layers + USE globalData,only:iname_watLayer ! named variable defining the total water state variable for snow+soil layers + USE globalData,only:iname_liqLayer ! named variable defining the liquid water state variable for snow+soil layers + USE globalData,only:iname_matLayer ! named variable defining the matric head state variable for soil layers + USE globalData,only:iname_lmpLayer ! named variable defining the liquid matric potential state variable for soil layers + + ! constants + USE multiconst,only:& + LH_fus, & ! latent heat of fusion (J kg-1) + iden_ice, & ! intrinsic density of ice (kg m-3) + iden_water ! intrinsic density of liquid water (kg m-3) + ! privacy + implicit none + private::printResidDAE + public::computResidSundials + contains + + ! ********************************************************************************************************** + ! public subroutine computResidSundials: compute the residual vector + ! ********************************************************************************************************** + subroutine computResidSundials(& + ! input: model control + nSnow, & ! intent(in): number of snow layers + nSoil, & ! intent(in): number of soil layers + nLayers, & ! intent(in): total number of layers + ! input: flux vectors + sMul, & ! intent(in): state vector multiplier (used in the residual calculations) + fVec, & ! intent(in): flux vector + ! input: state variables (already disaggregated into scalars and vectors) + scalarCanopyTempTrial, & ! intent(in): + mLayerTempTrial, & ! intent(in) + scalarCanairTempPrime, & ! intent(in): trial value for the temperature of the canopy air space (K) + scalarCanopyTempPrime, & ! intent(in): trial value for the temperature of the vegetation canopy (K) + scalarCanopyWatPrime, & ! + mLayerTempPrime, & ! intent(in): trial value for the temperature of each snow and soil layer (K) water content + scalarAquiferStoragePrime, & ! intent(in): trial value of storage of water in the aquifer (m) + ! input: diagnostic variables defining the liquid water and ice content (function of state variables) + scalarCanopyIcePrime, & ! intent(in): trial value for the ice on the vegetation canopy (kg m-2) + scalarCanopyLiqPrime, & ! intent(in): + mLayerVolFracIcePrime, & ! intent(in): trial value for the volumetric ice in each snow and soil layer (-) + mLayerVolFracWatPrime, & + mLayerVolFracLiqPrime, & + scalarCanopyCmTrial, & + mLayerCmTrial, & ! intent(in) + ! input: data structures + prog_data, & ! intent(in): model prognostic variables for a local HRU + diag_data, & ! intent(in): model diagnostic variables for a local HRU + flux_data, & ! intent(in): model fluxes for a local HRU + indx_data, & ! intent(in): index data + ! output + rAdd, & ! intent(out): additional (sink) terms on the RHS of the state equation + rVec, & ! intent(out): residual vector + err,message) ! intent(out): error control + ! -------------------------------------------------------------------------------------------------------------------------------- + implicit none + ! input: model control + integer(i4b),intent(in) :: nSnow ! number of snow layers + integer(i4b),intent(in) :: nSoil ! number of soil layers + integer(i4b),intent(in) :: nLayers ! total number of layers in the snow+soil domain + ! input: flux vectors + real(qp),intent(in) :: sMul(:) ! NOTE: qp ! state vector multiplier (used in the residual calculations) + real(rkind),intent(in) :: fVec(:) ! flux vector + ! input: state variables (already disaggregated into scalars and vectors) + real(rkind),intent(in) :: scalarCanopyTempTrial + real(rkind),intent(in) :: mLayerTempTrial(:) + real(rkind),intent(in) :: scalarCanairTempPrime ! trial value for temperature of the canopy air space (K) + real(rkind),intent(in) :: scalarCanopyTempPrime ! trial value for temperature of the vegetation canopy (K) + real(rkind),intent(in) :: scalarCanopyWatPrime ! derivative value for liquid water storage in the canopy (kg m-2) + real(rkind),intent(in) :: mLayerTempPrime(:) ! trial value for temperature of each snow/soil layer (K) content + real(rkind),intent(in) :: scalarAquiferStoragePrime ! trial value of aquifer storage (m) + ! input: diagnostic variables defining the liquid water and ice content (function of state variables) + real(rkind),intent(in) :: scalarCanopyIcePrime ! trial value for mass of ice on the vegetation canopy (kg m-2) + real(rkind),intent(in) :: scalarCanopyLiqPrime + real(rkind),intent(in) :: mLayerVolFracIcePrime(:) ! trial value for volumetric fraction of ice (-) + real(rkind),intent(in) :: mLayerVolFracLiqPrime(:) + real(rkind),intent(in) :: mLayerVolFracWatPrime(:) + real(qp),intent(in) :: scalarCanopyCmTrial + real(qp),intent(in) :: mLayerCmTrial(:) + ! input: data structures + type(var_dlength),intent(in) :: prog_data ! prognostic variables for a local HRU + type(var_dlength),intent(in) :: diag_data ! diagnostic variables for a local HRU + type(var_dlength),intent(in) :: flux_data ! model fluxes for a local HRU + type(var_ilength),intent(in) :: indx_data ! indices defining model states and layers + ! output + real(rkind),intent(out) :: rAdd(:) ! additional (sink) terms on the RHS of the state equation + real(qp),intent(out) :: rVec(:) ! NOTE: qp ! residual vector + integer(i4b),intent(out) :: err ! error code + character(*),intent(out) :: message ! error message + ! -------------------------------------------------------------------------------------------------------------------------------- + ! local variables + ! -------------------------------------------------------------------------------------------------------------------------------- + integer(i4b) :: iLayer ! index of layer within the snow+soil domain + integer(i4b),parameter :: ixVegVolume=1 ! index of the desired vegetation control volumne (currently only one veg layer) + real(rkind),dimension(nLayers) :: mLayerVolFracHydPrime ! vector of volumetric water content (-), either liquid water content or total water content + real(rkind) :: scalarCanopyHydPrime ! trial value for canopy water (kg m-2), either liquid water content or total water content + ! -------------------------------------------------------------------------------------------------------------------------------- + ! -------------------------------------------------------------------------------------------------------------------------------- + ! link to the necessary variables for the residual computations + associate(& + ! canopy and layer depth + canopyDepth => diag_data%var(iLookDIAG%scalarCanopyDepth)%dat(1) ,& ! intent(in): [dp] canopy depth (m) + mLayerDepth => prog_data%var(iLookPROG%mLayerDepth)%dat ,& ! intent(in): [dp(:)] depth of each layer in the snow-soil sub-domain (m) + ! model fluxes (sink terms in the soil domain) + mLayerTranspire => flux_data%var(iLookFLUX%mLayerTranspire)%dat ,& ! intent(in): [dp] transpiration loss from each soil layer (m s-1) + mLayerBaseflow => flux_data%var(iLookFLUX%mLayerBaseflow)%dat ,& ! intent(in): [dp(:)] baseflow from each soil layer (m s-1) + mLayerCompress => diag_data%var(iLookDIAG%mLayerCompress)%dat ,& ! intent(in): [dp(:)] change in storage associated with compression of the soil matrix (-) + ! number of state variables of a specific type + nSnowSoilNrg => indx_data%var(iLookINDEX%nSnowSoilNrg )%dat(1) ,& ! intent(in): [i4b] number of energy state variables in the snow+soil domain + nSnowSoilHyd => indx_data%var(iLookINDEX%nSnowSoilHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the snow+soil domain + nSoilOnlyHyd => indx_data%var(iLookINDEX%nSoilOnlyHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the soil domain + ! model indices + ixCasNrg => indx_data%var(iLookINDEX%ixCasNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy air space energy state variable + ixVegNrg => indx_data%var(iLookINDEX%ixVegNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy energy state variable + ixVegHyd => indx_data%var(iLookINDEX%ixVegHyd)%dat(1) ,& ! intent(in): [i4b] index of canopy hydrology state variable (mass) + ixAqWat => indx_data%var(iLookINDEX%ixAqWat)%dat(1) ,& ! intent(in): [i4b] index of water storage in the aquifer + ixSnowSoilNrg => indx_data%var(iLookINDEX%ixSnowSoilNrg)%dat ,& ! intent(in): [i4b(:)] indices for energy states in the snow+soil subdomain + ixSnowSoilHyd => indx_data%var(iLookINDEX%ixSnowSoilHyd)%dat ,& ! intent(in): [i4b(:)] indices for hydrology states in the snow+soil subdomain + ixSoilOnlyHyd => indx_data%var(iLookINDEX%ixSoilOnlyHyd)%dat ,& ! intent(in): [i4b(:)] indices for hydrology states in the soil subdomain + ixStateType => indx_data%var(iLookINDEX%ixStateType)%dat ,& ! intent(in): [i4b(:)] indices defining the type of the state (iname_nrgLayer...) + ixHydCanopy => indx_data%var(iLookINDEX%ixHydCanopy)%dat ,& ! intent(in): [i4b(:)] index of the hydrology states in the canopy domain + ixHydType => indx_data%var(iLookINDEX%ixHydType)%dat ,& ! intent(in): [i4b(:)] named variables defining the type of hydrology states in snow+soil domain + layerType => indx_data%var(iLookINDEX%layerType)%dat & ! intent(in): [i4b(:)] named variables defining the type of layer in snow+soil domain + ) ! association to necessary variables for the residual computations + ! -------------------------------------------------------------------------------------------------------------------------------- + ! initialize error control + err=0; message="computResidSundials/" + + ! --- + ! * compute sink terms... + ! ----------------------- + + ! intialize additional terms on the RHS as zero + rAdd(:) = 0._rkind + + ! compute energy associated with melt freeze for the vegetation canopy + if(ixVegNrg/=integerMissing) rAdd(ixVegNrg) = rAdd(ixVegNrg) + LH_fus*scalarCanopyIcePrime/canopyDepth ! energy associated with melt/freeze (J m-3) + ! compute energy associated with melt/freeze for snow + ! NOTE: allow expansion of ice during melt-freeze for snow; deny expansion of ice during melt-freeze for soil + if(nSnowSoilNrg>0)then + do concurrent (iLayer=1:nLayers,ixSnowSoilNrg(iLayer)/=integerMissing) ! (loop through non-missing energy state variables in the snow+soil domain) + select case( layerType(iLayer) ) + case(iname_snow); rAdd( ixSnowSoilNrg(iLayer) ) = rAdd( ixSnowSoilNrg(iLayer) ) + LH_fus*iden_ice * mLayerVolFracIcePrime(iLayer) + case(iname_soil); rAdd( ixSnowSoilNrg(iLayer) ) = rAdd( ixSnowSoilNrg(iLayer) ) + LH_fus*iden_water * mLayerVolFracIcePrime(iLayer) + end select + end do ! looping through non-missing energy state variables in the snow+soil domain + endif + + ! sink terms soil hydrology (-) + ! NOTE 1: state variable is volumetric water content, so melt-freeze is not included + ! NOTE 2: ground evaporation was already included in the flux at the upper boundary + ! NOTE 3: rAdd(ixSnowOnlyWat)=0, and is defined in the initialization above + ! NOTE 4: same sink terms for matric head and liquid matric potential + if(nSoilOnlyHyd>0)then + do concurrent (iLayer=1:nSoil,ixSoilOnlyHyd(iLayer)/=integerMissing) ! (loop through non-missing hydrology state variables in the snow+soil domain) + rAdd( ixSoilOnlyHyd(iLayer) ) = rAdd( ixSoilOnlyHyd(iLayer) ) + (mLayerTranspire(iLayer) - mLayerBaseflow(iLayer) )/mLayerDepth(iLayer+nSnow) - mLayerCompress(iLayer) + end do ! looping through non-missing energy state variables in the snow+soil domain + endif + + ! --- + ! * compute the residual vector... + ! -------------------------------- + + ! compute the residual vector for the vegetation canopy + ! NOTE: sMul(ixVegHyd) = 1, but include as it converts all variables to quadruple precision + ! --> energy balance + if(ixCasNrg/=integerMissing) rVec(ixCasNrg) = sMul(ixCasNrg)*scalarCanairTempPrime - ( fVec(ixCasNrg) + rAdd(ixCasNrg) ) + if(ixVegNrg/=integerMissing) rVec(ixVegNrg) = sMul(ixVegNrg) * scalarCanopyTempPrime + scalarCanopyCmTrial * scalarCanopyWatPrime/canopyDepth - ( fVec(ixVegNrg) + rAdd(ixVegNrg) ) + ! --> mass balance + if(ixVegHyd/=integerMissing)then + scalarCanopyHydPrime = merge(scalarCanopyWatPrime, scalarCanopyLiqPrime, (ixStateType( ixHydCanopy(ixVegVolume) )==iname_watCanopy) ) + rVec(ixVegHyd) = sMul(ixVegHyd)*scalarCanopyHydPrime - ( fVec(ixVegHyd) + rAdd(ixVegHyd) ) + endif + + ! compute the residual vector for the snow and soil sub-domains for energy + if(nSnowSoilNrg>0)then + do concurrent (iLayer=1:nLayers,ixSnowSoilNrg(iLayer)/=integerMissing) ! (loop through non-missing energy state variables in the snow+soil domain) + rVec( ixSnowSoilNrg(iLayer) ) = sMul( ixSnowSoilNrg(iLayer) ) * mLayerTempPrime(iLayer) + mLayerCmTrial(iLayer) * mLayerVolFracWatPrime(iLayer) - ( fVec( ixSnowSoilNrg(iLayer) ) + rAdd( ixSnowSoilNrg(iLayer) ) ) + end do ! looping through non-missing energy state variables in the snow+soil domain + endif + + ! compute the residual vector for the snow and soil sub-domains for hydrology + ! NOTE: residual depends on choice of state variable + if(nSnowSoilHyd>0)then + do concurrent (iLayer=1:nLayers,ixSnowSoilHyd(iLayer)/=integerMissing) ! (loop through non-missing hydrology state variables in the snow+soil domain) + ! (get the correct state variable) + mLayerVolFracHydPrime(iLayer) = merge(mLayerVolFracWatPrime(iLayer), mLayerVolFracLiqPrime(iLayer), (ixHydType(iLayer)==iname_watLayer .or. ixHydType(iLayer)==iname_matLayer) ) + ! (compute the residual) + rVec( ixSnowSoilHyd(iLayer) ) = mLayerVolFracHydPrime(iLayer) - ( fVec( ixSnowSoilHyd(iLayer) ) + rAdd( ixSnowSoilHyd(iLayer) ) ) + end do ! looping through non-missing energy state variables in the snow+soil domain + endif + + ! compute the residual vector for the aquifer + if(ixAqWat/=integerMissing) rVec(ixAqWat) = sMul(ixAqWat)*scalarAquiferStoragePrime - ( fVec(ixAqWat) + rAdd(ixAqWat) ) + + ! print result + if(globalPrintFlag)then + write(*,'(a,1x,100(e12.5,1x))') 'rVec = ', rVec(min(iJac1,size(rVec)):min(iJac2,size(rVec))) + write(*,'(a,1x,100(e12.5,1x))') 'fVec = ', fVec(min(iJac1,size(rVec)):min(iJac2,size(rVec))) + !print*, 'PAUSE:'; read(*,*) + endif + + !call printResidDAE(nSnow,nSoil,nLayers,indx_data,rAdd,rVec) + + ! check + if(any(isNan(rVec)))then + call printResidDAE(nSnow,nSoil,nLayers,indx_data,rAdd,rVec) + message=trim(message)//'we found NaN' + err=20; return + endif + + ! end association with the necessary variabiles for the residual calculations + end associate + + end subroutine computResidSundials + + ! ********************************************************************************************************** + ! private subroutine printResidDAE: print the residual vector mainly for debugging + ! ********************************************************************************************************** + subroutine printResidDAE( & + ! input: model control + nSnow, & ! intent(in): number of snow layers + nSoil, & ! intent(in): number of soil layers + nLayers, & ! intent(in): total number of layers + ! input: data structures + indx_data, & ! intent(in): index data + ! output + rAdd, & ! intent(out): additional (sink) terms on the RHS of the state equation + rVec) ! intent(out): residual vector + + ! -------------------------------------------------------------------------------------------------------------------------------- + implicit none + ! input: model control + integer(i4b),intent(in) :: nSnow ! number of snow layers + integer(i4b),intent(in) :: nSoil ! number of soil layers + integer(i4b),intent(in) :: nLayers ! total number of layers in the snow+soil domain + type(var_ilength),intent(in) :: indx_data ! indices defining model states and layers + ! output + real(rkind),intent(in) :: rAdd(:) ! additional (sink) terms on the RHS of the state equation + real(qp),intent(in) :: rVec(:) ! NOTE: qp ! residual vector + ! -------------------------------------------------------------------------------------------------------------------------------- + ! local variables + ! -------------------------------------------------------------------------------------------------------------------------------- + integer(i4b) :: iLayer ! index of layer within the snow+soil domain + ! -------------------------------------------------------------------------------------------------------------------------------- + ! link to the necessary variables for the residual computations + associate(& + ! number of state variables of a specific type + nSnowSoilNrg => indx_data%var(iLookINDEX%nSnowSoilNrg )%dat(1) ,& ! intent(in): [i4b] number of energy state variables in the snow+soil domain + nSnowSoilHyd => indx_data%var(iLookINDEX%nSnowSoilHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the snow+soil domain + nSoilOnlyHyd => indx_data%var(iLookINDEX%nSoilOnlyHyd )%dat(1) ,& ! intent(in): [i4b] number of hydrology variables in the soil domain + ! model indices + ixCasNrg => indx_data%var(iLookINDEX%ixCasNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy air space energy state variable + ixVegNrg => indx_data%var(iLookINDEX%ixVegNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy energy state variable + ixVegHyd => indx_data%var(iLookINDEX%ixVegHyd)%dat(1) ,& ! intent(in): [i4b] index of canopy hydrology state variable (mass) + ixAqWat => indx_data%var(iLookINDEX%ixAqWat)%dat(1) ,& ! intent(in): [i4b] index of water storage in the aquifer + ixSnowSoilNrg => indx_data%var(iLookINDEX%ixSnowSoilNrg)%dat ,& ! intent(in): [i4b(:)] indices for energy states in the snow+soil subdomain + ixSnowSoilHyd => indx_data%var(iLookINDEX%ixSnowSoilHyd)%dat ,& ! intent(in): [i4b(:)] indices for hydrology states in the snow+soil subdomain + ixSoilOnlyHyd => indx_data%var(iLookINDEX%ixSoilOnlyHyd)%dat ,& ! intent(in): [i4b(:)] indices for hydrology states in the soil subdomain + ixStateType => indx_data%var(iLookINDEX%ixStateType)%dat ,& ! intent(in): [i4b(:)] indices defining the type of the state (iname_nrgLayer...) + ixHydCanopy => indx_data%var(iLookINDEX%ixHydCanopy)%dat ,& ! intent(in): [i4b(:)] index of the hydrology states in the canopy domain + ixHydType => indx_data%var(iLookINDEX%ixHydType)%dat ,& ! intent(in): [i4b(:)] named variables defining the type of hydrology states in snow+soil domain + layerType => indx_data%var(iLookINDEX%layerType)%dat & ! intent(in): [i4b(:)] named variables defining the type of layer in snow+soil domain + ) ! association to necessary variables for the residual computations + ! -------------------------------------------------------------------------------------------------------------------------------- + + if(ixVegNrg/=integerMissing) print *, 'rAdd(ixVegNrg) = ', rAdd(ixVegNrg) + + if(nSnowSoilNrg>0)then + do concurrent (iLayer=1:nLayers,ixSnowSoilNrg(iLayer)/=integerMissing) + select case( layerType(iLayer) ) + case(iname_snow) + print *, 'rAdd( ixSnowSoilNrg(iLayer) ) = ', rAdd( ixSnowSoilNrg(iLayer) ) + case(iname_soil); print *, 'rAdd( ixSnowSoilNrg(iLayer) ) = ', rAdd( ixSnowSoilNrg(iLayer) ) + end select + end do + endif + + if(nSoilOnlyHyd>0)then + do concurrent (iLayer=1:nSoil,ixSoilOnlyHyd(iLayer)/=integerMissing) + print *, 'rAdd( ixSoilOnlyHyd(iLayer) ) = ', rAdd( ixSoilOnlyHyd(iLayer) ) + end do + endif + + if(ixCasNrg/=integerMissing) print *, 'rVec(ixCasNrg) = ', rVec(ixCasNrg) + if(ixVegNrg/=integerMissing) print *, 'rVec(ixVegNrg) = ', rVec(ixVegNrg) + if(ixVegHyd/=integerMissing)then + print *, 'rVec(ixVegHyd) = ', rVec(ixVegHyd) + endif + + if(nSnowSoilNrg>0)then + do concurrent (iLayer=1:nLayers,ixSnowSoilNrg(iLayer)/=integerMissing) + print *, 'rVec( ixSnowSoilNrg(iLayer) ) = ', rVec( ixSnowSoilNrg(iLayer) ) + end do + endif + + if(nSnowSoilHyd>0)then + do concurrent (iLayer=1:nLayers,ixSnowSoilHyd(iLayer)/=integerMissing) + print *, 'rVec( ixSnowSoilHyd(iLayer) ) = ', rVec( ixSnowSoilHyd(iLayer) ) + end do + endif + + if(ixAqWat/=integerMissing) print *, ' rVec(ixAqWat) = ', rVec(ixAqWat) + + end associate + + end subroutine printResidDAE + + end module computResidSundials_module + \ No newline at end of file diff --git a/utils/laugh_tests/colbeck1976/actors_out.txt b/utils/laugh_tests/colbeck1976/actors_out.txt new file mode 100644 index 0000000..3fe835e --- /dev/null +++ b/utils/laugh_tests/colbeck1976/actors_out.txt @@ -0,0 +1,349 @@ +SETTINGS FOR SUMMA_ACTOR +Output Structure Size = 250 +Max GRUs Per Job = 250 + +SETTINGS FOR JOB_ACTOR +File Manager Path = /Summa-Actors/utils/laugh_tests/colbeck1976/settings/summa_fileManager_colbeck1976-exp1.txt +output_csv = false +Job Actor Initalized + +----------File_Access_Actor Started---------- +Initalizing Output Structure +GRU Actor Has Started + +SETTINGS FOR HRU_ACTOR +Print Output = true +Print Output every 1 timesteps + +All Forcing Files Loaded +1 - Timestep = 1 +1 - Timestep = 2 +1 - Timestep = 3 +1 - Timestep = 4 +1 - Timestep = 5 +1 - Timestep = 6 +1 - Timestep = 7 +1 - Timestep = 8 +1 - Timestep = 9 +1 - Timestep = 10 +1 - Timestep = 11 + /Summa-Actors/utils/laugh_tests/colbeck1976/settings/summa_defineModelOutput.txt + HERE + +INFO: aspect not found in the input attribute file, continuing ... + + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 0.0000000000000000 + iden_water = 1000.0000000000000 + newSWE = 300.60000000000002 + delSWE = 0.60000000000002274 + massBalance = 2.2759572004815709E-014 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 4.9587017904807871E-017 + iden_water = 1000.0000000000000 + newSWE = 301.19999999999783 + delSWE = 0.59999999999780584 + massBalance = 7.8104189782379763E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 3.9772892257868010E-016 + iden_water = 1000.0000000000000 + newSWE = 301.79999999997375 + delSWE = 0.59999999997592113 + massBalance = -2.1516122217235534E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 1.3458248922576075E-015 + iden_water = 1000.0000000000000 + newSWE = 302.39999999989305 + delSWE = 0.59999999991930508 + massBalance = 5.4511950509095186E-014 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 3.1983853384144842E-015 + iden_water = 1000.0000000000000 + newSWE = 302.99999999970117 + delSWE = 0.59999999980811936 + massBalance = 2.2537527399890678E-014 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 6.2630410898092297E-015 + iden_water = 1000.0000000000000 + newSWE = 303.59999999932541 + delSWE = 0.59999999962423090 + massBalance = 1.3322676295501878E-014 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 1.0850552546761924E-014 + iden_water = 1000.0000000000000 + newSWE = 304.19999999867446 + delSWE = 0.59999999934905190 + massBalance = 8.5043083686286991E-014 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 1.7274818811323277E-014 + iden_water = 1000.0000000000000 + newSWE = 304.79999999763805 + delSWE = 0.59999999896359668 + massBalance = 8.5820239803524601E-014 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 2.5852886556788660E-014 + iden_water = 1000.0000000000000 + newSWE = 305.39999999608659 + delSWE = 0.59999999844853846 + massBalance = -2.8843594179761567E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 3.6904959605757223E-014 + iden_water = 1000.0000000000000 + newSWE = 305.99999999387285 + delSWE = 0.59999999778625579 + massBalance = 5.5333515547317802E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 5.0754385602336477E-014 + iden_water = 1000.0000000000000 + newSWE = 306.59999999082720 + delSWE = 0.59999999695435235 + massBalance = -3.8458125573015423E-013 +1 - Timestep = 12 +1 - Timestep = 13 +1 - Timestep = 14 +1 - Timestep = 15 +1 - Timestep = 16 +1 - Timestep = 17 +1 - Timestep = 18 +1 - Timestep = 19 +1 - Timestep = 20 +1 - Timestep = 21 +1 - Timestep = 22 +1 - Timestep = 23 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 6.7727690522882357E-014 + iden_water = 1000.0000000000000 + newSWE = 307.19999998676326 + delSWE = 0.59999999593605935 + massBalance = -2.7933211299568939E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 8.8154616957988339E-014 + iden_water = 1000.0000000000000 + newSWE = 307.79999998147423 + delSWE = 0.59999999471096999 + massBalance = 2.4702462297909733E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 1.1236808114754007E-013 + iden_water = 1000.0000000000000 + newSWE = 308.39999997473200 + delSWE = 0.59999999325776798 + massBalance = -1.4710455076283324E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 1.4070416063307484E-013 + iden_water = 1000.0000000000000 + newSWE = 308.99999996628952 + delSWE = 0.59999999155752448 + massBalance = -2.2593038551121936E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 1.7350215945493872E-013 + iden_water = 1000.0000000000000 + newSWE = 309.59999995587907 + delSWE = 0.59999998958954848 + massBalance = -3.2196467714129540E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 2.1110465794128700E-013 + iden_water = 1000.0000000000000 + newSWE = 310.19999994321324 + delSWE = 0.59999998733417215 + massBalance = 4.5163872641751368E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 2.5385736837767799E-013 + iden_water = 1000.0000000000000 + newSWE = 310.79999992798213 + delSWE = 0.59999998476888550 + massBalance = 3.2762681456688370E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 3.0210926319812266E-013 + iden_water = 1000.0000000000000 + newSWE = 311.39999990985564 + delSWE = 0.59999998187350911 + massBalance = 6.4837024638109142E-014 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 3.5621257403774553E-013 + iden_water = 1000.0000000000000 + newSWE = 311.99999988848282 + delSWE = 0.59999997862718146 + massBalance = -6.4170890823334048E-014 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 4.1652276488339124E-013 + iden_water = 1000.0000000000000 + newSWE = 312.59999986349135 + delSWE = 0.59999997500852942 + massBalance = -1.0469403122215226E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 4.8339858993693743E-013 + iden_water = 1000.0000000000000 + newSWE = 313.19999983448764 + delSWE = 0.59999997099629354 + massBalance = 2.0894397323445446E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 5.5720191270652321E-013 + iden_water = 1000.0000000000000 + newSWE = 313.79999980105532 + delSWE = 0.59999996656767962 + massBalance = -2.0561330416057899E-013 +1 - Timestep = 24 +1 - Timestep = 25 +1 - Timestep = 26 +1 - Timestep = 27 +1 - Timestep = 28 +1 - Timestep = 29 +1 - Timestep = 30 +1 - Timestep = 31 +1 - Timestep = 32 +1 - Timestep = 33 +1 - Timestep = 34 +1 - Timestep = 35 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 6.3829798560848913E-013 + iden_water = 1000.0000000000000 + newSWE = 314.39999976275755 + delSWE = 0.59999996170222403 + massBalance = 1.0313971898767704E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 7.2705534112198951E-013 + iden_water = 1000.0000000000000 + newSWE = 314.99999971913411 + delSWE = 0.59999995637656411 + massBalance = -1.1546319456101628E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 8.2384578888474908E-013 + iden_water = 1000.0000000000000 + newSWE = 315.59999966970321 + delSWE = 0.59999995056909938 + massBalance = -1.5332179970073412E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 9.2904443946227717E-013 + iden_water = 1000.0000000000000 + newSWE = 316.19999961396070 + delSWE = 0.59999994425749037 + massBalance = 1.5665246877460959E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 1.0430297137078835E-012 + iden_water = 1000.0000000000000 + newSWE = 316.79999955137873 + delSWE = 0.59999993741803337 + massBalance = -1.8385293287792592E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 1.1661833323718264E-012 + iden_water = 1000.0000000000000 + newSWE = 317.39999948140763 + delSWE = 0.59999993002890051 + massBalance = -9.9587005308876542E-014 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 1.2988901893932967E-012 + iden_water = 1000.0000000000000 + newSWE = 317.99999940347448 + delSWE = 0.59999992206684283 + massBalance = 2.5424107263916085E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 1.4415385941891714E-012 + iden_water = 1000.0000000000000 + newSWE = 318.59999931698201 + delSWE = 0.59999991350753135 + massBalance = -1.5309975509580909E-013 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 1.5945202341461499E-012 + iden_water = 1000.0000000000000 + newSWE = 319.19999922131086 + delSWE = 0.59999990432885397 + massBalance = 6.8056671409522096E-014 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 1.7582301232143246E-012 + iden_water = 1000.0000000000000 + newSWE = 319.79999911581712 + delSWE = 0.59999989450625435 + massBalance = 6.1728400169158704E-014 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 1.9330666721022761E-012 + iden_water = 1000.0000000000000 + newSWE = 320.39999899983303 + delSWE = 0.59999988401591509 + massBalance = -8.4598994476436928E-014 + effSnowfall = 0.0000000000000000 + effRainfall = 1.0000000000000000E-002 + averageSnowSublimation = 0.0000000000000000 + averageSnowDrainage = 2.1194312945133812E-012 + iden_water = 1000.0000000000000 + newSWE = 320.99999887266722 + delSWE = 0.59999987283418932 + massBalance = 6.6946448384896939E-014 +1 - Timestep = 36 +1 - Timestep = 37 +1 - Timestep = 38 +1 - Timestep = 39 +1 - Timestep = 40 +1 - Timestep = 41 +1 - Timestep = 42 +1 - Timestep = 43 +1 - Timestep = 44 +Error: RunPhysics - HRU = 1 - indxGRU = 1 - refGRU = 1 - Timestep = 44 +*** unexpected message [id: 10, name: user.scheduled-actor]: message(run_failure(), 11@38DDCD1CBB03D559767271D57A7B827F2487A98D#36110, 1, 20) -- GitLab