diff --git a/build/source/engine/sundials/eval8DAE.f90 b/build/source/engine/sundials/eval8DAE.f90 deleted file mode 100644 index 6ed2414b97d1a0068aa6db9b1a87d389e7e3615f..0000000000000000000000000000000000000000 --- a/build/source/engine/sundials/eval8DAE.f90 +++ /dev/null @@ -1,757 +0,0 @@ - -module eval8DAE_module - -! data types -USE nrtype - -! access missing values -USE globalData,only:integerMissing ! missing integer -USE globalData,only:realMissing ! missing double precision number -USE globalData,only:quadMissing ! missing quadruple precision number - -! access the global print flag -USE globalData,only:globalPrintFlag - -! 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 -USE globalData,only:model_decisions ! model decision structure - - -! 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) - -! provide access to the 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) - zlookup, & - model_options ! defines the model decisions - -! indices that define elements of the data structures -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 - -! 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 - enthalpyFD ! heat capacity using enthalpy - -! 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 - -implicit none -private -public::eval8DAE - -contains - -! ********************************************************************************************************** -! public subroutine eval8DAE: compute the residual vector -! ********************************************************************************************************** -subroutine eval8DAE(& - ! input: model control - dt_cur, & ! intent(in): current stepsize - dt, & ! intent(in): entire time step - nSnow, & ! intent(in): number of snow layers - nSoil, & ! intent(in): number of soil layers - nLayers, & ! intent(in): total number of layers - nState, & ! intent(in): total number of state variables - checkFeas, & ! intent(in): flag to indicate if we are checking for feasibility - firstSubStep, & ! intent(in): flag to indicate if we are processing the first sub-step - firstFluxCall, & ! intent(inout) flag to indicate if we are processing the first flux call - firstSplitOper, & ! intent(inout) flag to indicate if we are processing the first flux call in a splitting operation - computeVegFlux, & ! intent(in): flag to indicate if we need to compute fluxes over vegetation - scalarSolution, & ! intent(in): flag to indicate the scalar solution - requireLWBal, & ! intent(in): flag to indicate if we need longwave to be balanced - ! input: state vectors - stateVec, & ! intent(in): model state vector - stateVecPrime, & ! intent(in): derivative of model state vector - sMul, & ! intent(inout): state vector multiplier (used in the residual calculations) - ! input: data structures - model_decisions, & ! intent(in): model decisions - lookup_data, & ! intent(in): lookup data - type_data, & ! intent(in): type of vegetation and soil - attr_data, & ! intent(in): spatial attributes - mpar_data, & ! intent(in): model parameters - forc_data, & ! intent(in): model forcing data - bvar_data, & ! intent(in): average model variables for the entire basin - 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 - flux_data, & ! intent(inout): model fluxes for a local HRU - deriv_data, & ! intent(inout): derivatives in model fluxes w.r.t. relevant state variables - ! input-output: - dBaseflow_dMatric, & ! intent(out): derivative in baseflow w.r.t. matric head (s-1) - scalarCanopyTempTrial, & ! intent(out): trial value of canopy temperature (K) - scalarCanopyTempPrev, & ! intent(in): value of canopy temperature (K) - scalarCanopyIceTrial, & ! intent(out): trial value for mass of ice on the vegetation canopy (kg m-2) - scalarCanopyIcePrev, & ! intent(in): value for mass of ice on the vegetation canopy (kg m-2) - scalarCanopyLiqTrial, & ! intent(out): trial value of canopy liquid water (kg m-2) - scalarCanopyLiqPrev, & ! intent(in): value of canopy liquid water (kg m-2) - scalarCanopyEnthalpyTrial,& ! intent(out): trial value for enthalpy of the vegetation canopy (J m-3) - scalarCanopyEnthalpyPrev,& ! intent(in): value for enthalpy of the vegetation canopy (J m-3) - mLayerTempTrial, & ! intent(out): trial vector of layer temperature (K) - mLayerTempPrev, & ! intent(in): vector of layer temperature (K) - mLayerMatricHeadLiqTrial,& ! intent(out): trial value for liquid water matric potential (m) - mLayerMatricHeadTrial, & ! intent(out): trial value for total water matric potential (m) - mLayerMatricHeadPrev, & ! intent(in): value for total water matric potential (m) - mLayerVolFracWatTrial, & ! intent(out): trial vector of volumetric total water content (-) - mLayerVolFracWatPrev, & ! intent(in): vector of volumetric total water content (-) - mLayerVolFracIceTrial, & ! intent(out): trial vector of volumetric ice water content (-) - mLayerVolFracIcePrev, & ! intent(in): vector of volumetric ice water content (-) - mLayerVolFracLiqTrial, & ! intent(out): trial vector of volumetric liquid water content (-) - mLayerVolFracLiqPrev, & ! intent(in): vector of volumetric liquid water content (-) - scalarAquiferStorageTrial,& ! intent(out): trial value of storage of water in the aquifer (m) - scalarAquiferStoragePrev,& ! intent(in): value of storage of water in the aquifer (m) - mLayerEnthalpyPrev, & ! intent(in): vector of enthalpy for snow+soil layers (J m-3) - mLayerEnthalpyTrial, & ! intent(out): trial vector of enthalpy for snow+soil layers (J m-3) - ixSaturation, & ! intent(inout): index of the lowest saturated layer - feasible, & ! intent(out): flag to denote the feasibility of the solution - fluxVec, & ! intent(out): flux vector - resSink, & ! intent(out): additional (sink) terms on the RHS of the state equation - resVec, & ! intent(out): residual vector - err,message) ! intent(out): error control - ! -------------------------------------------------------------------------------------------------------------------------------- - ! provide access to subroutines - USE varExtrSundials_module, only:varExtract2 ! extract variables from the state vector - USE varExtrSundials_module, only:varExtractSundials - USE updateVarsSundials_module, only:updateVarsSundials ! update variables - USE t2enthalpy_module, only:t2enthalpy_T ! compute enthalpy - USE computFlux_module, only:soilCmpresSundials ! compute soil compression - USE computFlux_module, only:computFlux ! compute fluxes given a state vector - USE computHeatCap_module,only:computHeatCap ! compute heat capacity - USE computHeatCap_module,only:computHeatCapAnalytic ! compute heat capacity - USE computHeatCap_module,only:computCm - USE computHeatCap_module, only:computStatMult - USE computResidDAE_module,only:computResidDAE ! compute residuals given a state vector - USE computThermConduct_module,only:computThermConduct - USE computEnthalpy_module,only:computEnthalpy - USE computEnthalpy_module,only:computEnthalpyPrime - implicit none - ! -------------------------------------------------------------------------------------------------------------------------------- - ! -------------------------------------------------------------------------------------------------------------------------------- - ! input: model control - real(rkind),intent(in) :: dt_cur - 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(kind=8),intent(in) :: nState ! total number of state variables - logical(lgt),intent(in) :: checkFeas ! flag to indicate if we are checking for feasibility - logical(lgt),intent(in) :: firstSubStep ! flag to indicate if we are processing the first sub-step - logical(lgt),intent(inout) :: firstFluxCall - logical(lgt),intent(inout) :: firstSplitOper ! flag to indicate if we are processing the first flux call in a splitting operation - 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 - logical(lgt),intent(in) :: requireLWBal ! flag to indicate if we need longwave to be balanced - ! input: state vectors - real(rkind),intent(in) :: stateVec(:) ! model state vector - real(rkind),intent(in) :: stateVecPrime(:) ! model state vector - real(qp),intent(inout) :: sMul(:) ! NOTE: qp ! state vector multiplier (used in the residual calculations) - ! input: data structures - type(model_options),intent(in) :: model_decisions(:) ! model decisions - type(zLookup),intent(in) :: lookup_data ! lookup tables - type(var_i), intent(in) :: type_data ! type of vegetation and soil - type(var_d), intent(in) :: attr_data ! spatial attributes - type(var_dlength), intent(in) :: mpar_data ! model parameters - type(var_d), intent(in) :: forc_data ! model forcing data - type(var_dlength), intent(in) :: bvar_data ! model variables for the local basin - 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) :: flux_data ! model fluxes 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(out) :: dBaseflow_dMatric(:,:) ! derivative in baseflow w.r.t. matric head (s-1) - ! output: flux and residual vectors - real(rkind),intent(out) :: scalarCanopyTempTrial ! trial value for temperature of the vegetation canopy (K) - real(rkind),intent(in) :: scalarCanopyTempPrev ! previous value for temperature of the vegetation canopy (K) - real(rkind),intent(out) :: scalarCanopyIceTrial ! trial value for mass of ice on the vegetation canopy (kg m-2) - real(rkind),intent(in) :: scalarCanopyIcePrev ! previous value for mass of ice on the vegetation canopy (kg m-2) - real(rkind),intent(out) :: scalarCanopyLiqTrial ! trial value for mass of ice on the vegetation canopy (kg m-2) - real(rkind),intent(in) :: scalarCanopyLiqPrev ! previous value for mass of ice on the vegetation canopy (kg m-2) - real(rkind),intent(out) :: scalarCanopyEnthalpyTrial ! trial value for enthalpy of the vegetation canopy (J m-3) - real(rkind),intent(in) :: scalarCanopyEnthalpyPrev ! previous value of enthalpy of the vegetation canopy (J m-3) - real(rkind),intent(out) :: mLayerTempTrial(:) ! trial vector of layer temperature (K) - real(rkind),intent(in) :: mLayerTempPrev(:) - real(rkind),intent(out) :: mLayerMatricHeadLiqTrial(:) ! trial value for liquid water matric potential (m) - real(rkind),intent(out) :: mLayerMatricHeadTrial(:) ! trial value for total water matric potential (m) - real(rkind),intent(in) :: mLayerMatricHeadPrev(:) ! value for total water matric potential (m) - real(rkind),intent(out) :: mLayerVolFracWatTrial(:) ! trial vector of volumetric total water content (-) - real(rkind),intent(in) :: mLayerVolFracWatPrev(:) ! vector of volumetric total water content (-) - real(rkind),intent(out) :: mLayerVolFracIceTrial(:) ! trial vector of volumetric ice water content (-) - real(rkind),intent(in) :: mLayerVolFracIcePrev(:) ! vector of volumetric ice water content (-) - real(rkind),intent(out) :: mLayerVolFracLiqTrial(:) ! trial vector of volumetric liquid water content (-) - real(rkind),intent(in) :: mLayerVolFracLiqPrev(:) ! vector of volumetric liquid water content (-) - real(rkind),intent(out) :: scalarAquiferStorageTrial ! trial value of storage of water in the aquifer (m) - real(rkind),intent(in) :: scalarAquiferStoragePrev ! value of storage of water in the aquifer (m) - real(rkind),intent(in) :: mLayerEnthalpyPrev(:) ! vector of enthalpy for snow+soil layers (J m-3) - real(rkind),intent(out) :: mLayerEnthalpyTrial(:) ! trial vector of enthalpy for snow+soil layers (J m-3) - integer(i4b),intent(inout) :: ixSaturation ! index of the lowest saturated layer - logical(lgt),intent(out) :: feasible ! flag to denote the feasibility of the solution - real(rkind),intent(out) :: fluxVec(:) ! flux vector - real(rkind),intent(out) :: resSink(:) ! sink terms on the RHS of the flux equation - real(qp),intent(out) :: resVec(:) ! NOTE: qp ! residual vector - ! 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) :: scalarCanopyWatTrial ! trial value for liquid water storage in the canopy (kg m-2) - ! 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 (-) - ! enthalpy - real(rkind) :: scalarCanairEnthalpy ! enthalpy of the canopy air space (J m-3) - ! other local variables - integer(i4b) :: iLayer ! index of model layer in the snow+soil domain - integer(i4b) :: jState(1) ! index of model state for the scalar solution within the soil domain - integer(i4b) :: ixBeg,ixEnd ! index of indices for the soil compression routine - integer(i4b),parameter :: ixVegVolume=1 ! index of the desired vegetation control volumne (currently only one veg layer) - real(rkind) :: xMin,xMax ! minimum and maximum values for water content - real(rkind),parameter :: canopyTempMax=500._rkind ! expected maximum value for the canopy temperature (K) - character(LEN=256) :: cmessage ! error message of downwind routine - real(rkind) :: scalarCanopyCmTrial ! trial value of Cm for the canopy - real(rkind),dimension(nLayers) :: mLayerCmTrial ! trial vector of Cm for snow+soil - logical(lgt),parameter :: updateCp=.true. ! flag to indicate if we update Cp at each step - logical(lgt),parameter :: needCm=.false. ! flag to indicate if the energy equation contains Cm = dH_T/dTheta_m - - - - ! -------------------------------------------------------------------------------------------------------------------------------- - ! 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 - ! snow parameters - snowfrz_scale => mpar_data%var(iLookPARAM%snowfrz_scale)%dat(1) ,& ! intent(in): [dp] scaling parameter for the snow freezing curve (K-1) - ! 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) - theta_res => mpar_data%var(iLookPARAM%theta_res)%dat ,& ! intent(in): [dp(:)] residual volumetric water content (-) - ! 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 state variables - scalarSfcMeltPond => prog_data%var(iLookPROG%scalarSfcMeltPond)%dat(1) ,& ! intent(in): [dp] ponded water caused by melt of the "snow without a layer" (kg m-2) - mLayerVolFracIce => prog_data%var(iLookPROG%mLayerVolFracIce)%dat ,& ! intent(in): [dp(:)] volumetric fraction of ice (-) - ! soil compression - scalarSoilCompress => diag_data%var(iLookDIAG%scalarSoilCompress)%dat(1) ,& ! intent(in): [dp] total change in storage associated with compression of the soil matrix (kg m-2) - mLayerCompress => diag_data%var(iLookDIAG%mLayerCompress)%dat ,& ! intent(in): [dp(:)] change in storage associated with compression of the soil matrix (-) - ! derivatives - dVolTot_dPsi0 => deriv_data%var(iLookDERIV%dVolTot_dPsi0)%dat ,& ! intent(in): [dp(:)] 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 (m-1) - ! mapping - ixMapFull2Subset => indx_data%var(iLookINDEX%ixMapFull2Subset)%dat ,& ! intent(in): [i4b(:)] mapping of full state vector to the state subset - ixControlVolume => indx_data%var(iLookINDEX%ixControlVolume)%dat ,& ! intent(in): [i4b(:)] index of control volume for different domains (veg, snow, soil) - ! indices - ixCasNrg => indx_data%var(iLookINDEX%ixCasNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy air space energy state variable (nrg) - ixVegNrg => indx_data%var(iLookINDEX%ixVegNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy energy state variable (nrg) - ixVegHyd => indx_data%var(iLookINDEX%ixVegHyd)%dat(1) ,& ! intent(in): [i4b] index of canopy hydrology state variable (mass) - ixSnowOnlyNrg => indx_data%var(iLookINDEX%ixSnowOnlyNrg)%dat ,& ! intent(in): [i4b(:)] indices for energy states in the snow subdomain - ixSnowSoilHyd => indx_data%var(iLookINDEX%ixSnowSoilHyd)%dat ,& ! intent(in): [i4b(:)] indices for hydrology states in the snow+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(:)] index of the type of hydrology states in snow+soil domain - layerType => indx_data%var(iLookINDEX%layerType)%dat ,& ! intent(in): [i4b(:)] layer type (iname_soil or iname_snow) - heatCapVegTrial => diag_data%var(iLookDIAG%scalarBulkVolHeatCapVeg)%dat(1) ,& ! intent(out): volumetric heat capacity of vegetation canopy - mLayerHeatCapTrial => diag_data%var(iLookDIAG%mLayerVolHtCapBulk)%dat & ! intent(out): heat capacity for snow and soil - ) ! association to variables in the data structures - ! -------------------------------------------------------------------------------------------------------------------------------- - ! initialize error control - err=0; message="eval8DAE/" - feasible=.true. - - ! check the feasibility of the solution - if (checkFeas) then - ! check that the canopy air space temperature is reasonable - if(ixCasNrg/=integerMissing)then - if(stateVec(ixCasNrg) > canopyTempMax) feasible=.false. - if(stateVec(ixCasNrg) > canopyTempMax) message=trim(message)//'canopy air space temp high,' - if(.not.feasible) write(*,'(a,1x,L1,1x,10(f20.10,1x))') 'feasible, max, stateVec( ixCasNrg )', feasible, canopyTempMax, stateVec(ixCasNrg) - endif - - ! check that the canopy temperature is reasonable - if(ixVegNrg/=integerMissing)then - if(stateVec(ixVegNrg) > canopyTempMax) feasible=.false. - if(stateVec(ixVegNrg) > canopyTempMax) message=trim(message)//'canopy temp high,' - if(.not.feasible) write(*,'(a,1x,L1,1x,10(f20.10,1x))') 'feasible, max, stateVec( ixVegNrg )', feasible, canopyTempMax, stateVec(ixVegNrg) - endif - - ! check canopy liquid water is not negative - if(ixVegHyd/=integerMissing)then - if(stateVec(ixVegHyd) < 0._rkind) feasible=.false. - if(stateVec(ixVegHyd) < 0._rkind) message=trim(message)//'canopy water negative,' - if(.not.feasible) write(*,'(a,1x,L1,1x,10(f20.10,1x))') 'feasible, min, stateVec( ixVegHyd )', feasible, 0._rkind, stateVec(ixVegHyd) - end if - - ! check snow temperature is below freezing - if(count(ixSnowOnlyNrg/=integerMissing)>0)then - if(any(stateVec( pack(ixSnowOnlyNrg,ixSnowOnlyNrg/=integerMissing) ) > Tfreeze)) feasible=.false. - if(any(stateVec( pack(ixSnowOnlyNrg,ixSnowOnlyNrg/=integerMissing) ) > Tfreeze)) message=trim(message)//'snow temp above freezing,' - do iLayer=1,nSnow - if(.not.feasible) write(*,'(a,1x,i4,1x,L1,1x,10(f20.10,1x))') 'iLayer, feasible, max, stateVec( ixSnowOnlyNrg(iLayer) )', iLayer, feasible, Tfreeze, stateVec( ixSnowOnlyNrg(iLayer) ) - end do - endif - - ! loop through non-missing hydrology state variables in the snow+soil domain - do concurrent (iLayer=1:nLayers,ixSnowSoilHyd(iLayer)/=integerMissing) - - ! check the minimum and maximum water constraints - if(ixHydType(iLayer)==iname_watLayer .or. ixHydType(iLayer)==iname_liqLayer)then - - ! --> minimum - if (layerType(iLayer) == iname_soil) then - xMin = theta_res(iLayer-nSnow) - else - xMin = 0._rkind - endif - - ! --> maximum - select case( layerType(iLayer) ) - case(iname_snow); xMax = merge(iden_ice, 1._rkind - mLayerVolFracIce(iLayer), ixHydType(iLayer)==iname_watLayer) - case(iname_soil); xMax = merge(theta_sat(iLayer-nSnow), theta_sat(iLayer-nSnow) - mLayerVolFracIce(iLayer), ixHydType(iLayer)==iname_watLayer) - end select - - ! --> check - if(stateVec( ixSnowSoilHyd(iLayer) ) < xMin .or. stateVec( ixSnowSoilHyd(iLayer) ) > xMax) feasible=.false. - if(stateVec( ixSnowSoilHyd(iLayer) ) < xMin .or. stateVec( ixSnowSoilHyd(iLayer) ) > xMax) message=trim(message)//'layer water outside bounds,' - if(.not.feasible) write(*,'(a,1x,i4,1x,L1,1x,10(f20.10,1x))') 'iLayer, feasible, stateVec( ixSnowSoilHyd(iLayer) ), xMin, xMax = ', iLayer, feasible, stateVec( ixSnowSoilHyd(iLayer) ), xMin, xMax - - endif ! if water states - - end do ! loop through non-missing hydrology state variables in the snow+soil domain - - ! early return for non-feasible solutions - if(.not.feasible)then - fluxVec(:) = realMissing - resVec(:) = quadMissing - message=trim(message)//'non-feasible' - err=20; return - end if - - end if ! ( feasibility check ) - - ! get the start and end indices for the soil compression calculations - if(scalarSolution)then - jState = pack(ixControlVolume, ixMapFull2Subset/=integerMissing) - ixBeg = jState(1) - ixEnd = jState(1) - else - ixBeg = 1 - ixEnd = nSoil - endif - - ! initialize to state variable from the last update - scalarCanopyTempTrial = scalarCanopyTempPrev - scalarCanopyLiqTrial = scalarCanopyLiqPrev - scalarCanopyIceTrial = scalarCanopyIcePrev - mLayerTempTrial = mLayerTempPrev - mLayerVolFracWatTrial = mLayerVolFracWatPrev - mLayerVolFracLiqTrial = mLayerVolFracLiqPrev - mLayerVolFracIceTrial = mLayerVolFracIcePrev - mLayerMatricHeadTrial = mLayerMatricHeadPrev - scalarAquiferStorageTrial = scalarAquiferStoragePrev - - ! extract variables from the model state vector - call varExtract2(& - ! input - stateVec, & ! intent(in): model state vector (mixed units) - indx_data, & ! intent(in): indices defining model states and layers - ! output: variables for the vegetation canopy - scalarCanairTempTrial, & ! intent(out): trial value of canopy air temperature (K) - scalarCanopyTempTrial, & ! intent(out): trial value of canopy temperature (K) - scalarCanopyWatTrial, & ! intent(out): trial value of canopy total water (kg m-2) - scalarCanopyLiqTrial, & ! intent(out): trial value of canopy liquid water (kg m-2) - ! output: variables for the snow-soil domain - mLayerTempTrial, & ! intent(out): trial vector of layer temperature (K) - mLayerVolFracWatTrial, & ! intent(out): trial vector of volumetric total water content (-) - mLayerVolFracLiqTrial, & ! intent(out): trial vector of volumetric liquid water content (-) - mLayerMatricHeadTrial, & ! intent(out): trial vector of total water matric potential (m) - mLayerMatricHeadLiqTrial, & ! intent(out): trial vector of liquid water matric potential (m) - ! output: variables for the aquifer - scalarAquiferStorageTrial,& ! intent(out): trial value 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 varExtractSundials(& - ! input - stateVecPrime, & ! intent(in): derivative of model state vector (mixed units) - indx_data, & ! intent(in): indices defining model states and layers - ! output: variables for the vegetation canopy - scalarCanairTempPrime, & ! intent(out): derivative of canopy air temperature (K) - scalarCanopyTempPrime, & ! intent(out): derivative of canopy temperature (K) - scalarCanopyWatPrime, & ! intent(out): derivative of canopy total water (kg m-2) - scalarCanopyLiqPrime, & ! intent(out): derivative of canopy liquid water (kg m-2) - ! output: variables for the snow-soil domain - mLayerTempPrime, & ! intent(out): derivative of layer temperature (K) - mLayerVolFracWatPrime, & ! intent(out): derivative of volumetric total water content (-) - mLayerVolFracLiqPrime, & ! intent(out): derivative of volumetric liquid water content (-) - mLayerMatricHeadPrime, & ! intent(out): derivative of total water matric potential (m) - mLayerMatricHeadLiqPrime, & ! intent(out): derivative of liquid water matric potential (m) - ! output: variables for the aquifer - scalarAquiferStoragePrime,& ! intent(out): 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_cur, & - .false., & ! intent(in): logical flag to adjust temperature to account for the energy - 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 - mLayerVolFracWatPrev, & ! intent(in) - mLayerMatricHeadPrev, & ! intent(in) - 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) - - - ! print the water content - if(globalPrintFlag)then - if(iJac1<nSnow) write(*,'(a,10(f16.10,1x))') 'mLayerVolFracWatTrial = ', mLayerVolFracWatTrial(iJac1:min(iJac2,nSnow)) - if(iJac1<nSnow) write(*,'(a,10(f16.10,1x))') 'mLayerVolFracLiqTrial = ', mLayerVolFracLiqTrial(iJac1:min(iJac2,nSnow)) - if(iJac1<nSnow) write(*,'(a,10(f16.10,1x))') 'mLayerVolFracIceTrial = ', mLayerVolFracIceTrial(iJac1:min(iJac2,nSnow)) - endif - - - if(updateCp)then - ! *** compute volumetric heat capacity C_p - if(model_decisions(iLookDECISIONS%howHeatCap)%iDecision == enthalpyFD)then - ! compute H_T - call t2enthalpy_T(& - ! input: data structures - diag_data, & ! intent(in): model diagnostic variables for a local HRU - mpar_data, & ! intent(in): parameter data structure - indx_data, & ! intent(in): model indices - lookup_data, & ! intent(in): lookup table data structure - ! input: state variables for the vegetation canopy - scalarCanairTempTrial, & ! intent(in): trial value of canopy air temperature (K) - scalarCanopyTempTrial, & ! intent(in): trial value of canopy temperature (K) - scalarCanopyWatTrial, & ! intent(in): trial value of canopy total water (kg m-2) - scalarCanopyIceTrial, & ! intent(in): trial value for canopy ice content (kg m-2) - ! input: variables for the snow-soil domain - mLayerTempTrial, & ! intent(in): trial vector of layer temperature (K) - mLayerVolFracWatTrial, & ! intent(in): trial vector of volumetric total water content (-) - mLayerMatricHeadTrial, & ! intent(in): trial vector of total water matric potential (m) - mLayerVolFracIceTrial, & ! intent(in): trial vector of volumetric fraction of ice (-) - ! output: enthalpy - scalarCanairEnthalpy, & ! intent(out): enthalpy of the canopy air space (J m-3) - scalarCanopyEnthalpyTrial, & ! intent(out): enthalpy of the vegetation canopy (J m-3) - mLayerEnthalpyTrial, & ! intent(out): enthalpy of each snow+soil layer (J m-3) - ! output: error control - err,cmessage) ! intent(out): error control - if(err/=0)then; message=trim(message)//trim(cmessage); return; endif - - ! *** compute volumetric heat capacity C_p = dH_T/dT - call computHeatCap(& - ! input: control variables - nLayers, & ! intent(in): number of layers (soil+snow) - computeVegFlux, & ! intent(in): flag to denote if computing the vegetation flux - canopyDepth, & ! intent(in): canopy depth (m) - ! input data structures - mpar_data, & ! intent(in): model parameters - indx_data, & ! intent(in): model layer indices - ! input: state variables - scalarCanopyIceTrial, & ! intent(in): trial value for canopy ice content (kg m-2) - scalarCanopyLiqTrial, & ! intent(in): trial value for the liquid water on the vegetation canopy (kg m-2) - scalarCanopyTempTrial, & ! intent(in): trial value of canopy temperature (K) - scalarCanopyTempPrev, & ! intent(in): previous value of canopy temperature (K) - scalarCanopyEnthalpyTrial, & ! intent(in): trial enthalpy of the vegetation canopy (J m-3) - scalarCanopyEnthalpyPrev, & ! intent(in): previous enthalpy of the vegetation canopy (J m-3) - mLayerVolFracIceTrial, & ! intent(in): volumetric fraction of ice at the start of the sub-step (-) - mLayerVolFracLiqTrial, & ! intent(in): volumetric fraction of liquid water at the start of the sub-step (-) - mLayerTempTrial, & ! intent(in): trial temperature - mLayerTempPrev, & ! intent(in): previous temperature - mLayerEnthalpyTrial, & ! intent(in): trial enthalpy for snow and soil - mLayerEnthalpyPrev, & ! intent(in): previous enthalpy for snow and soil - ! output - heatCapVegTrial, & ! intent(out): volumetric heat capacity of vegetation canopy - mLayerHeatCapTrial, & ! intent(out): heat capacity for snow and soil - ! output: error control - err,cmessage) ! intent(out): error control - if(err/=0)then; message=trim(message)//trim(cmessage); return; endif - ! to conserve energy compute finite difference approximation of (theta_ice)' - if(dt_cur > 1e-14_rkind) then - scalarCanopyIcePrime = ( scalarCanopyIceTrial - scalarCanopyIcePrev ) / dt_cur - do concurrent (iLayer=1:nLayers) - mLayerVolFracIcePrime(iLayer) = ( mLayerVolFracIceTrial(iLayer) - mLayerVolFracIcePrev(iLayer) ) / dt_cur - end do - endif ! if dt_cur is not too samll - else ! if using closed formula of heat capacity - call computHeatCapAnalytic(& - ! input: control variables - computeVegFlux, & ! intent(in): flag to denote if computing the vegetation flux - canopyDepth, & ! intent(in): canopy depth (m) - ! input: state variables - scalarCanopyIceTrial, & ! intent(in) - scalarCanopyLiqTrial, & ! intent(in) - mLayerVolFracIceTrial, & ! intent(in): volumetric fraction of ice at the start of the sub-step (-) - mLayerVolFracLiqTrial, & ! intent(in): fraction of liquid water at the start of the sub-step (-) - ! input data structures - mpar_data, & ! intent(in): model parameters - indx_data, & ! intent(in): model layer indices - ! output - heatCapVegTrial, & ! intent(out): volumetric heat capacity of vegetation canopy - mLayerHeatCapTrial, & ! intent(out): volumetric heat capacity of soil and snow - ! output: error control - err,cmessage) ! intent(out): error control - endif - - ! compute multiplier of state vector - call computStatMult(& - ! input - heatCapVegTrial, & ! intent(in): volumetric heat capacity of vegetation canopy - mLayerHeatCapTrial, & ! intent(in): volumetric heat capacity of soil and snow - diag_data, & ! intent(in): model diagnostic variables for a local HRU - indx_data, & ! intent(in): indices defining model states and layers - ! output - sMul, & ! intent(out): multiplier for state vector (used in the residual calculations) - err,cmessage) ! intent(out): error control - if(err/=0)then; message=trim(message)//trim(cmessage); return; endif ! (check for errors) - - ! update thermal conductivity - call computThermConduct(& - ! input: control variables - computeVegFlux, & ! intent(in): flag to denote if computing the vegetation flux - ! input: state variables - scalarCanopyIceTrial, & ! intent(in) - scalarCanopyLiqTrial, & ! intent(in) - mLayerVolFracIceTrial, & ! intent(in): volumetric fraction of ice at the start of the sub-step (-) - mLayerVolFracLiqTrial, & ! intent(in): volumetric fraction of liquid water at the start of the sub-step (-) - ! input/output: data structures - mpar_data, & ! intent(in): model parameters - indx_data, & ! intent(in): model layer indices - prog_data, & ! intent(in): model prognostic variables for a local HRU - diag_data, & ! intent(inout): model diagnostic variables for a local HRU - err,cmessage) ! intent(out): error control - if(err/=0)then; err=55; message=trim(message)//trim(cmessage); return; end if - - end if ! updateCp - - - if(needCm)then - ! compute C_m - call computCm(& - ! input: control variables - computeVegFlux, & ! intent(in): flag to denote if computing the vegetation flux - ! input: state variables - scalarCanopyTempTrial, & ! intent(in) - mLayerTempTrial, & ! intent(in): volumetric fraction of liquid water at the start of the sub-step (-) - mLayerMatricHeadTrial, & ! intent(in) - ! input data structures - mpar_data, & ! intent(in): model parameters - indx_data, & ! intent(in): model layer indices - ! output - scalarCanopyCmTrial, & ! intent(out): Cm for vegetation - mLayerCmTrial, & ! intent(out): Cm for soil and snow - err,cmessage) ! intent(out): error control - else - scalarCanopyCmTrial = 0._qp - mLayerCmTrial = 0._qp - end if ! needCm - - - ! save the number of flux calls per time step - indx_data%var(iLookINDEX%numberFluxCalc)%dat(1) = indx_data%var(iLookINDEX%numberFluxCalc)%dat(1) + 1 - ! compute the fluxes for a given state vector - call computFlux(& - ! input-output: model control - nSnow, & ! intent(in): number of snow layers - nSoil, & ! intent(in): number of soil layers - nLayers, & ! intent(in): total number of layers - firstSubStep, & ! intent(in): flag to indicate if we are processing the first sub-step - firstFluxCall, & ! intent(inout): flag to denote the first flux call - firstSplitOper, & ! intent(in): flag to indicate if we are processing the first flux call in a splitting operation - computeVegFlux, & ! intent(in): flag to indicate if we need to compute fluxes over vegetation - scalarSolution, & ! intent(in): flag to indicate the scalar solution - requireLWBal, & ! intent(in): flag to indicate if we need longwave to be balanced - scalarSfcMeltPond/dt, & ! intent(in): drainage from the surface melt pond (kg m-2 s-1) - ! input: state variables - scalarCanairTempTrial, & ! intent(in): trial value for the temperature of the canopy air space (K) - scalarCanopyTempTrial, & ! intent(in): trial value for the temperature of the vegetation canopy (K) - mLayerTempTrial, & ! intent(in): trial value for the temperature of each snow and soil layer (K) - mLayerMatricHeadLiqTrial, & ! intent(in): trial value for the liquid water matric potential in each soil layer (m) - mLayerMatricHeadTrial, & ! intent(in): trial vector of total water matric potential (m) - scalarAquiferStorageTrial, & ! intent(in): trial value of storage of water in the aquifer (m) - ! input: diagnostic variables defining the liquid water and ice content - scalarCanopyLiqTrial, & ! intent(in): trial value for the liquid water on the vegetation canopy (kg m-2) - scalarCanopyIceTrial, & ! intent(in): trial value for the ice on the vegetation canopy (kg m-2) - mLayerVolFracLiqTrial, & ! intent(in): trial value for the volumetric liquid water content in each snow and soil layer (-) - mLayerVolFracIceTrial, & ! intent(in): trial value for the volumetric ice in each snow and soil layer (-) - ! input: data structures - model_decisions, & ! intent(in): model decisions - type_data, & ! intent(in): type of vegetation and soil - attr_data, & ! intent(in): spatial attributes - mpar_data, & ! intent(in): model parameters - forc_data, & ! intent(in): model forcing data - bvar_data, & ! intent(in): average model variables for the entire basin - prog_data, & ! intent(in): model prognostic variables for a local HRU - indx_data, & ! intent(in): index data - ! input-output: data structures - diag_data, & ! intent(inout): model diagnostic variables for a local HRU - flux_data, & ! intent(inout): model fluxes for a local HRU - deriv_data, & ! intent(out): derivatives in model fluxes w.r.t. relevant state variables - ! input-output: flux vector and baseflow derivatives - ixSaturation, & ! intent(inout): index of the lowest saturated layer (NOTE: only computed on the first iteration) - dBaseflow_dMatric, & ! intent(out): derivative in baseflow w.r.t. matric head (s-1), we will use it later in computeJacobSundials - fluxVec, & ! intent(out): flux vector (mixed units) - ! 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) - - firstSplitOper = .true. - - - ! compute soil compressibility (-) and its derivative w.r.t. matric head (m) - ! NOTE: we already extracted trial matrix head and volumetric liquid water as part of the flux calculations - call soilCmpresSundials(& - ! input: - ixRichards, & ! intent(in): choice of option for Richards' equation - ixBeg,ixEnd, & ! intent(in): start and end indices defining desired layers - mLayerMatricHeadPrime(1:nSoil), & ! intent(in): matric head at the start of the time step (m) - mLayerVolFracLiqTrial(nSnow+1:nLayers), & ! intent(in): trial value for the volumetric liquid water content in each soil layer (-) - mLayerVolFracIceTrial(nSnow+1:nLayers), & ! intent(in): trial value for the volumetric ice content in each soil layer (-) - specificStorage, & ! intent(in): specific storage coefficient (m-1) - theta_sat, & ! intent(in): soil porosity (-) - ! output: - mLayerCompress, & ! intent(inout): compressibility of the soil matrix (-) - dCompress_dPsi, & ! intent(inout): derivative in compressibility w.r.t. matric head (m-1) - err,cmessage) ! intent(out): error code and error message - if(err/=0)then; message=trim(message)//trim(cmessage); return; end if ! (check for errors) - - ! print *, 'dt = ', dt - ! print *, 'dt_cur = ', dt_cur - - ! compute the residual vector - call 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) - fluxVec, & ! intent(in): flux vector - ! input: state variables (already disaggregated into scalars and vectors) - scalarCanopyTempTrial, & ! intent(in): - mLayerTempTrial, & ! intent(in) - scalarCanairTempPrime, & ! intent(in): Prime value for the temperature of the canopy air space (K) - scalarCanopyTempPrime, & ! intent(in): Prime value for the temperature of the vegetation canopy (K) - scalarCanopyWatPrime, & - mLayerTempPrime, & ! intent(in): Prime value for the temperature of each snow and soil layer (K) - scalarAquiferStoragePrime, & ! intent(in): Prime 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): Prime value for the ice on the vegetation canopy (kg m-2) - scalarCanopyLiqPrime, & ! intent(in): - mLayerVolFracIcePrime, & ! intent(in): Prime value for the volumetric ice in each snow and soil layer (-) - mLayerVolFracWatPrime, & - mLayerVolFracLiqPrime, & - scalarCanopyCmTrial, & ! intent(in) Cm of vegetation canopy - mLayerCmTrial, & ! intent(in) Cm of soil and snow - ! 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 - resSink, & ! intent(out): additional (sink) terms on the RHS of the state equation - resVec, & ! intent(out): residual vector - err,cmessage) ! intent(out): error control - if(err/=0)then; message=trim(message)//trim(cmessage); return; end if ! (check for errors) - - !print *, '=====================================================================================' - - - ! end association with the information in the data structures - end associate - - -end subroutine eval8DAE -end module eval8DAE_module diff --git a/build/source/engine/sundials/evalDAE4IDA.f90 b/build/source/engine/sundials/evalDAE4IDA.f90 deleted file mode 100644 index 55b20b02ef2be29de45ff0e7326d43be6bf36ad0..0000000000000000000000000000000000000000 --- a/build/source/engine/sundials/evalDAE4IDA.f90 +++ /dev/null @@ -1,167 +0,0 @@ - - -module evalDAE4IDA_module - - - !======= Inclusions =========== - use, intrinsic :: iso_c_binding - use nrtype - use type4IDA - USE globalData,only:model_decisions ! model decision structure - USE globalData,only:flux_meta ! metadata on the model fluxes - ! provide access to the 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 - USE multiconst,only:iden_water ! intrinsic density of liquid water (kg m-3) - USE var_lookup,only:iLookDIAG - USE var_lookup,only:iLookPROG - USE var_lookup,only:iLookINDEX ! named variables for structure elements - USE var_lookup,only:iLookDERIV ! named variables for structure elements - - - ! privacy - implicit none - private - public::evalDAE4IDA - - -contains - - ! ********************************************************************************************************** - ! public function evalDAE4IDA: compute the residual vector F(t,y,y') required for IDA solver - ! ********************************************************************************************************** - ! Return values: - ! 0 = success, - ! 1 = recoverable error, - ! -1 = non-recoverable error - ! ---------------------------------------------------------------- - integer(c_int) function evalDAE4IDA(tres, sunvec_y, sunvec_yp, sunvec_r, user_data) & - result(ierr) bind(C,name='evalDAE4IDA') - - !======= Inclusions =========== - use, intrinsic :: iso_c_binding - use fida_mod - use fsundials_nvector_mod - use fnvector_serial_mod - use nrtype - use type4IDA - use eval8DAE_module,only:eval8DAE - - !======= Declarations ========= - implicit none - - ! calling variables - real(c_double), value :: tres ! current time t - type(N_Vector) :: sunvec_y ! solution N_Vector y - type(N_Vector) :: sunvec_yp ! derivative N_Vector y' - type(N_Vector) :: sunvec_r ! residual N_Vector F(t,y,y') - type(c_ptr), value :: user_data ! user-defined data - - - ! pointers to data in SUNDIALS vectors - type(eqnsData), pointer :: eqns_data ! equations data - real(rkind), pointer :: stateVec(:) - real(rkind), pointer :: stateVecPrime(:) - real(rkind), pointer :: rVec(:) - logical(lgt) :: feasible - integer(i4b) :: retval - real(c_double) :: stepsize_next(1) - !======= 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) - - retval = FIDAGetCurrentStep(eqns_data%ida_mem, stepsize_next) - if (retval /= 0) then - print *, 'Error in FIDAGetCurrentStep, retval = ', retval, '; halting' - stop 1 - end if - - ! compute the flux and the residual vector for a given state vector - call eval8DAE(& - ! input: model control - stepsize_next(1), & ! intent(in): current stepsize - 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%nState, & ! intent(in): number of state variables in the current subset - .false., & ! intent(in): do not check for feasibility inside Sundials loop - eqns_data%firstSubStep, & ! intent(in): flag to indicate if we are processing the first sub-step - eqns_data%firstFluxCall, & ! intent(inout): flag to indicate if we are processing the first flux call - eqns_data%firstSplitOper, & ! intent(inout): flag to indicate if we are processing the first flux call in a splitting operation - 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 - .false., & ! intent(in): do not require that longwave is balanced inside Sundials loop - ! input: state vectors - stateVec, & ! intent(in): model state vector - stateVecPrime, & ! intent(in): model state vector - eqns_data%sMul, & ! intent(inout): state vector multiplier (used in the residual calculations) - ! input: data structures - model_decisions, & ! intent(in): model decisions - eqns_data%lookup_data, & ! intent(in): lookup data - eqns_data%type_data, & ! intent(in): type of vegetation and soil - eqns_data%attr_data, & ! intent(in): spatial attributes - eqns_data%mpar_data, & ! intent(in): model parameters - eqns_data%forc_data, & ! intent(in): model forcing data - eqns_data%bvar_data, & ! intent(in): average model variables for the entire basin - 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%flux_data, & ! intent(inout): model fluxes for a local HRU (initial flux structure) - eqns_data%deriv_data, & ! intent(inout): derivatives in model fluxes w.r.t. relevant state variables - ! input-output: baseflow - eqns_data%dBaseflow_dMatric, & ! intent(out): derivative in baseflow w.r.t. matric head (s-1), we will use it later for Jacobian - eqns_data%scalarCanopyTempTrial, & ! intent(in): trial value of canopy temperature (K) - eqns_data%scalarCanopyTempPrev, & ! intent(in): previous value of canopy temperature (K) - eqns_data%scalarCanopyIceTrial, & ! intent(out): trial value for mass of ice on the vegetation canopy (kg m-2) - eqns_data%scalarCanopyIcePrev, & ! intent(in): value for mass of ice on the vegetation canopy (kg m-2) - eqns_data%scalarCanopyLiqTrial, & ! intent(out): trial value of canopy liquid water (kg m-2) - eqns_data%scalarCanopyLiqPrev, & ! intent(in): value of canopy liquid water (kg m-2) - eqns_data%scalarCanopyEnthalpyTrial,& ! intent(out): trial value for enthalpy of the vegetation canopy (J m-3) - eqns_data%scalarCanopyEnthalpyPrev, & ! intent(in): value for enthalpy of the vegetation canopy (J m-3) - eqns_data%mLayerTempTrial, & ! intent(out): trial vector of layer temperature (K) - eqns_data%mLayerTempPrev, & ! intent(in): vector of layer temperature (K) - eqns_data%mLayerMatricHeadLiqTrial,& ! intent(out): trial value for liquid water matric potential (m) - eqns_data%mLayerMatricHeadTrial, & ! intent(out): trial value for total water matric potential (m) - eqns_data%mLayerMatricHeadPrev, & ! intent(in): value for total water matric potential (m) - eqns_data%mLayerVolFracWatTrial, & ! intent(out): trial vector of volumetric total water content (-) - eqns_data%mLayerVolFracWatPrev, & ! intent(in): vector of volumetric total water content (-) - eqns_data%mLayerVolFracIceTrial, & ! intent(out): trial vector of volumetric ice water content (-) - eqns_data%mLayerVolFracIcePrev, & ! intent(in): vector of volumetric ice water content (-) - eqns_data%mLayerVolFracLiqTrial, & ! intent(out): trial vector of volumetric liquid water content (-) - eqns_data%mLayerVolFracLiqPrev, & ! intent(in): vector of volumetric liquid water content (-) - eqns_data%scalarAquiferStorageTrial, & ! intent(out): trial value of storage of water in the aquifer (m) - eqns_data%scalarAquiferStoragePrev, & ! intent(in): value of storage of water in the aquifer (m) - eqns_data%mLayerEnthalpyPrev, & ! intent(in): vector of enthalpy for snow+soil layers (J m-3) - eqns_data%mLayerEnthalpyTrial, & ! intent(out): trial vector of enthalpy for snow+soil layers (J m-3) - eqns_data%ixSaturation, & ! intent(inout): index of the lowest saturated layer - ! output - feasible, & ! intent(out): flag to denote the feasibility of the solution - eqns_data%fluxVec, & ! intent(out): flux vector - eqns_data%resSink, & ! intent(out): additional (sink) terms on the RHS of the state equation - rVec, & ! intent(out): residual vector - 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 - if(.not.feasible)then; eqns_data%message=trim(eqns_data%message)//'state vector not feasible'; ierr = 1; return; endif - - ! return success - ierr = 0 - return - - end function evalDAE4IDA - - -end module evalDAE4IDA_module diff --git a/build/source/engine/sundials/evalJac4IDA.f90 b/build/source/engine/sundials/evalJac4IDA.f90 deleted file mode 100644 index c7da698497ae19c7563c281f7e9a17a7e2fdcf94..0000000000000000000000000000000000000000 --- a/build/source/engine/sundials/evalJac4IDA.f90 +++ /dev/null @@ -1,128 +0,0 @@ - - -module evalJac4IDA_module - - - !======= Inclusions =========== - use, intrinsic :: iso_c_binding - use nrtype - use type4IDA - USE globalData,only:model_decisions ! model decision structure - ! provide access to the 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 - - - - ! privacy - implicit none - private - public::evalJac4IDA - - -contains - - ! ********************************************************************************************************** - ! public function evalJac4IDA: 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 evalJac4IDA(t, cj, sunvec_y, sunvec_yp, sunvec_r, & - sunmat_J, user_data, sunvec_temp1, sunvec_temp2, sunvec_temp3) & - result(ierr) bind(C,name='evalJac4IDA') - - !======= 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 - use eval8JacDAE_module,only:eval8JacDAE ! compute Jacobian matrix - !======= Declarations ========= - implicit none - - ! calling variables - real(c_double), value :: t ! current time - real(c_double), 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 eval8JacDAE(& - ! 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 evalJac4IDA - - -end module evalJac4IDA_module diff --git a/build/source/engine/sundials/updateVars4JacDAE.f90 b/build/source/engine/sundials/updateVars4JacDAE.f90 deleted file mode 100644 index 2f0d2e39248d0c28b2bcff8a4b2227b52f4e1673..0000000000000000000000000000000000000000 --- a/build/source/engine/sundials/updateVars4JacDAE.f90 +++ /dev/null @@ -1,816 +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 updateVars4JacDAE_module - -! data types -USE nrtype - -! missing values -USE globalData,only:integerMissing ! missing integer -USE globalData,only:realMissing ! missing real number - -! access the global print flag -USE globalData,only:globalPrintFlag - -! domain types -USE globalData,only:iname_cas ! named variables for canopy air space -USE globalData,only:iname_veg ! named variables for vegetation canopy -USE globalData,only:iname_snow ! named variables for snow -USE globalData,only:iname_soil ! named variables for soil -USE globalData,only:iname_aquifer ! named variables for the aquifer - -! 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 total water on the vegetation canopy -USE globalData,only:iname_liqCanopy ! named variable defining the mass of liquid 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 - -! metadata for information in the data structures -USE globalData,only:indx_meta ! metadata for the variables in the index structure - -! constants -USE multiconst,only:& - gravity, & ! acceleration of gravity (m s-2) - Tfreeze, & ! temperature at freezing (K) - Cp_air, & ! specific heat of air (J kg-1 K-1) - Cp_ice, & ! specific heat of ice (J kg-1 K-1) - Cp_water, & ! specific heat of liquid water (J kg-1 K-1) - LH_fus, & ! latent heat of fusion (J kg-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) - -! provide access to the 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) - -! provide access to indices that define elements of the data structures -USE var_lookup,only:iLookDIAG ! named variables for structure elements -USE var_lookup,only:iLookPROG ! named variables for structure elements -USE var_lookup,only:iLookDERIV ! named variables for structure elements -USE var_lookup,only:iLookPARAM ! named variables for structure elements -USE var_lookup,only:iLookINDEX ! named variables for structure elements - -! provide access to routines to update states -USE updatStateSundials_module,only:updateVegSundials ! update snow states -USE updatStateSundials_module,only:updateSnowSundials ! update snow states -USE updatStateSundials_module,only:updateSoilSundials2 ! update soil states - -! provide access to functions for the constitutive functions and derivatives -USE snow_utils_module,only:fracliquid ! compute the fraction of liquid water (snow) -USE snow_utils_module,only:dFracLiq_dTk ! differentiate the freezing curve w.r.t. temperature (snow) -USE soil_utils_module,only:dTheta_dTk ! differentiate the freezing curve w.r.t. temperature (soil) -USE soil_utils_module,only:dTheta_dPsi ! derivative in the soil water characteristic (soil) -USE soil_utils_module,only:dPsi_dTheta ! derivative in the soil water characteristic (soil) -USE soil_utils_module,only:matricHead ! compute the matric head based on volumetric water content -USE soil_utils_module,only:volFracLiq ! compute volumetric fraction of liquid water -USE soil_utils_module,only:crit_soilT ! compute critical temperature below which ice exists -USE soil_utilsSundials_module,only:liquidHeadSundials ! compute the liquid water matric potential -USE soil_utilsSundials_module,only:d2Theta_dPsi2 -USE soil_utilsSundials_module,only:d2Theta_dTk2 - -! IEEE checks -USE, intrinsic :: ieee_arithmetic ! check values (NaN, etc.) - -implicit none -private -public::updateVars4JacDAE - -contains - - ! ********************************************************************************************************** - ! public subroutine updateVars4JacDAE: compute diagnostic variables - ! ********************************************************************************************************** - subroutine updateVars4JacDAE(& - ! input - do_adjustTemp, & ! 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 - 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, & ! reza - mLayerVolFracWatPrime, & ! reza - mLayerVolFracLiqPrime, & ! reza - mLayerVolFracIcePrime, & ! reza - mLayerMatricHeadPrime, & ! reza - mLayerMatricHeadLiqPrime, & ! reza - ! output: error control - err,message) ! intent(out): error control - ! -------------------------------------------------------------------------------------------------------------------------------- - ! -------------------------------------------------------------------------------------------------------------------------------- - implicit none - ! input - logical(lgt) ,intent(in) :: do_adjustTemp ! flag to adjust temperature to account for the energy used in melt+freeze - type(var_dlength),intent(in) :: mpar_data ! model parameters for a local HRU - 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(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 - ! output: variables for the vegetation canopy - real(rkind),intent(inout) :: scalarCanopyTempTrial ! trial value of canopy temperature (K) - real(rkind),intent(inout) :: scalarCanopyWatTrial ! trial value of canopy total water (kg m-2) - real(rkind),intent(inout) :: scalarCanopyLiqTrial ! trial value of canopy liquid water (kg m-2) - real(rkind),intent(inout) :: scalarCanopyIceTrial ! trial value of canopy ice content (kg m-2) - real(rkind),intent(inout) :: scalarCanopyTempPrime ! trial value of canopy temperature (K) - real(rkind),intent(inout) :: scalarCanopyWatPrime ! trial value of canopy total water (kg m-2) - real(rkind),intent(inout) :: scalarCanopyLiqPrime ! trial value of canopy liquid water (kg m-2) - real(rkind),intent(inout) :: scalarCanopyIcePrime ! trial value of canopy ice content (kg m-2) - ! output: variables for the snow-soil domain - real(rkind),intent(inout) :: mLayerTempTrial(:) ! trial vector of layer temperature (K) - real(rkind),intent(inout) :: mLayerVolFracWatTrial(:) ! trial vector of volumetric total water content (-) - real(rkind),intent(inout) :: mLayerVolFracLiqTrial(:) ! trial vector of volumetric liquid water content (-) - real(rkind),intent(inout) :: mLayerVolFracIceTrial(:) ! trial vector of volumetric ice water content (-) - real(rkind),intent(inout) :: mLayerMatricHeadTrial(:) ! trial vector of total water matric potential (m) - real(rkind),intent(inout) :: mLayerMatricHeadLiqTrial(:) ! trial vector of liquid water matric potential (m) - - real(rkind),intent(inout) :: mLayerTempPrime(:) - real(rkind),intent(inout) :: mLayerVolFracWatPrime(:) ! reza - real(rkind),intent(inout) :: mLayerVolFracLiqPrime(:) ! reza - real(rkind),intent(inout) :: mLayerVolFracIcePrime(:) ! reza - real(rkind),intent(inout) :: mLayerMatricHeadPrime(:) ! reza - real(rkind),intent(inout) :: mLayerMatricHeadLiqPrime(:) ! reza - ! output: error control - integer(i4b),intent(out) :: err ! error code - character(*),intent(out) :: message ! error message - ! -------------------------------------------------------------------------------------------------------------------------------- - ! general local variables - integer(i4b) :: iState ! index of model state variable - integer(i4b) :: iLayer ! index of layer within the snow+soil domain - integer(i4b) :: ixFullVector ! index within full state vector - integer(i4b) :: ixDomainType ! name of a given model domain - integer(i4b) :: ixControlIndex ! index within a given model domain - integer(i4b) :: ixOther,ixOtherLocal ! index of the coupled state variable within the (full, local) vector - logical(lgt) :: isCoupled ! .true. if a given variable shared another state variable in the same control volume - logical(lgt) :: isNrgState ! .true. if a given variable is an energy state - logical(lgt),allocatable :: computedCoupling(:) ! .true. if computed the coupling for a given state variable - real(rkind) :: scalarVolFracLiq ! volumetric fraction of liquid water (-) - real(rkind) :: scalarVolFracIce ! volumetric fraction of ice (-) - real(rkind) :: scalarVolFracLiqPrime ! volumetric fraction of liquid water (-) - real(rkind) :: scalarVolFracIcePrime ! volumetric fraction of ice (-) - real(rkind) :: Tcrit ! critical soil temperature below which ice exists (K) - real(rkind) :: xTemp ! temporary temperature (K) - real(rkind) :: fLiq ! fraction of liquid water (-) - real(rkind) :: effSat ! effective saturation (-) - real(rkind) :: avPore ! available pore space (-) - character(len=256) :: cMessage ! error message of downwind routine - logical(lgt),parameter :: printFlag=.false. ! flag to turn on printing - ! iterative solution for temperature - real(rkind) :: meltNrg ! energy for melt+freeze (J m-3) - real(rkind) :: residual ! residual in the energy equation (J m-3) - real(rkind) :: derivative ! derivative in the energy equation (J m-3 K-1) - real(rkind) :: tempInc ! iteration increment (K) - integer(i4b) :: iter ! iteration index - integer(i4b) :: niter ! number of iterations - integer(i4b),parameter :: maxiter=100 ! maximum number of iterations - real(rkind),parameter :: nrgConvTol=1.e-4_rkind ! convergence tolerance for energy (J m-3) - real(rkind),parameter :: tempConvTol=1.e-6_rkind ! convergence tolerance for temperature (K) - real(rkind) :: critDiff ! temperature difference from critical (K) - real(rkind) :: tempMin ! minimum bracket for temperature (K) - real(rkind) :: tempMax ! maximum bracket for temperature (K) - logical(lgt) :: bFlag ! flag to denote that iteration increment was constrained using bi-section - real(rkind),parameter :: epsT=1.e-7_rkind ! small interval above/below critical temperature (K) - ! -------------------------------------------------------------------------------------------------------------------------------- - ! make association with variables in the data structures - associate(& - ! number of model layers, and layer type - nSnow => indx_data%var(iLookINDEX%nSnow)%dat(1) ,& ! intent(in): [i4b] total number of snow layers - nSoil => indx_data%var(iLookINDEX%nSoil)%dat(1) ,& ! intent(in): [i4b] total number of soil layers - nLayers => indx_data%var(iLookINDEX%nLayers)%dat(1) ,& ! intent(in): [i4b] total number of snow and soil layers - mLayerDepth => prog_data%var(iLookPROG%mLayerDepth)%dat ,& ! intent(in): [dp(:)] depth of each layer in the snow-soil sub-domain (m) - ! indices defining model states and layers - 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) - ! indices in the full vector for specific domains - ixNrgCanair => indx_data%var(iLookINDEX%ixNrgCanair)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in canopy air space domain - ixNrgCanopy => indx_data%var(iLookINDEX%ixNrgCanopy)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in the canopy domain - ixHydCanopy => indx_data%var(iLookINDEX%ixHydCanopy)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for hydrology states in the canopy domain - 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 - ! mapping between the full state vector and the state subset - ixMapFull2Subset => indx_data%var(iLookINDEX%ixMapFull2Subset)%dat ,& ! intent(in): [i4b(:)] list of indices in the state subset for each state in the full state vector - ixMapSubset2Full => indx_data%var(iLookINDEX%ixMapSubset2Full)%dat ,& ! intent(in): [i4b(:)] [state subset] list of indices of the full state vector in the state subset - ! type of domain, type of state variable, and index of control volume within domain - ixDomainType_subset => indx_data%var(iLookINDEX%ixDomainType_subset)%dat ,& ! intent(in): [i4b(:)] [state subset] id of domain for desired model state variables - ixControlVolume => indx_data%var(iLookINDEX%ixControlVolume)%dat ,& ! intent(in): [i4b(:)] index of the control volume for different domains (veg, snow, soil) - ixStateType => indx_data%var(iLookINDEX%ixStateType)%dat ,& ! intent(in): [i4b(:)] indices defining the type of the state (iname_nrgLayer...) - ! snow parameters - snowfrz_scale => mpar_data%var(iLookPARAM%snowfrz_scale)%dat(1) ,& ! intent(in): [dp] scaling parameter for the snow freezing curve (K-1) - ! depth-varying model parameters - vGn_m => diag_data%var(iLookDIAG%scalarVGn_m)%dat ,& ! intent(in): [dp(:)] van Genutchen "m" parameter (-) - vGn_n => mpar_data%var(iLookPARAM%vGn_n)%dat ,& ! intent(in): [dp(:)] van Genutchen "n" parameter (-) - vGn_alpha => mpar_data%var(iLookPARAM%vGn_alpha)%dat ,& ! intent(in): [dp(:)] van Genutchen "alpha" parameter (m-1) - theta_sat => mpar_data%var(iLookPARAM%theta_sat)%dat ,& ! intent(in): [dp(:)] soil porosity (-) - theta_res => mpar_data%var(iLookPARAM%theta_res)%dat ,& ! intent(in): [dp(:)] soil residual volumetric water content (-) - ! model diagnostic variables (heat capacity) - canopyDepth => diag_data%var(iLookDIAG%scalarCanopyDepth)%dat(1) ,& ! intent(in): [dp ] canopy depth (m) - scalarBulkVolHeatCapVeg => diag_data%var(iLookDIAG%scalarBulkVolHeatCapVeg)%dat(1),& ! intent(in): [dp ] volumetric heat capacity of the vegetation (J m-3 K-1) - mLayerVolHtCapBulk => diag_data%var(iLookDIAG%mLayerVolHtCapBulk)%dat ,& ! intent(in): [dp(:)] volumetric heat capacity in each layer (J m-3 K-1) - ! model diagnostic variables (fraction of liquid water) - scalarFracLiqVeg => diag_data%var(iLookDIAG%scalarFracLiqVeg)%dat(1) ,& ! intent(out): [dp] fraction of liquid water on vegetation (-) - mLayerFracLiqSnow => diag_data%var(iLookDIAG%mLayerFracLiqSnow)%dat ,& ! intent(out): [dp(:)] fraction of liquid water in each snow layer (-) - ! model states for the vegetation canopy - scalarCanairTemp => prog_data%var(iLookPROG%scalarCanairTemp)%dat(1) ,& ! intent(in): [dp] temperature of the canopy air space (K) - scalarCanopyTemp => prog_data%var(iLookPROG%scalarCanopyTemp)%dat(1) ,& ! intent(in): [dp] temperature of the vegetation canopy (K) - scalarCanopyWat => prog_data%var(iLookPROG%scalarCanopyWat)%dat(1) ,& ! intent(in): [dp] mass of total water on the vegetation canopy (kg m-2) - ! model state variable vectors for the snow-soil layers - mLayerTemp => prog_data%var(iLookPROG%mLayerTemp)%dat ,& ! intent(in): [dp(:)] temperature of each snow/soil layer (K) - mLayerVolFracWat => prog_data%var(iLookPROG%mLayerVolFracWat)%dat ,& ! intent(in): [dp(:)] volumetric fraction of total water (-) - mLayerMatricHead => prog_data%var(iLookPROG%mLayerMatricHead)%dat ,& ! intent(in): [dp(:)] total water matric potential (m) - mLayerMatricHeadLiq => diag_data%var(iLookDIAG%mLayerMatricHeadLiq)%dat ,& ! intent(in): [dp(:)] liquid water matric potential (m) - ! model diagnostic variables from a previous solution - scalarCanopyLiq => prog_data%var(iLookPROG%scalarCanopyLiq)%dat(1) ,& ! intent(in): [dp(:)] mass of liquid water on the vegetation canopy (kg m-2) - scalarCanopyIce => prog_data%var(iLookPROG%scalarCanopyIce)%dat(1) ,& ! intent(in): [dp(:)] mass of ice on the vegetation canopy (kg m-2) - mLayerVolFracLiq => prog_data%var(iLookPROG%mLayerVolFracLiq)%dat ,& ! intent(in): [dp(:)] volumetric fraction of liquid water (-) - mLayerVolFracIce => prog_data%var(iLookPROG%mLayerVolFracIce)%dat ,& ! intent(in): [dp(:)] volumetric fraction of ice (-) - ! derivatives - dVolTot_dPsi0 => deriv_data%var(iLookDERIV%dVolTot_dPsi0 )%dat ,& ! intent(out): [dp(:)] derivative in total water content w.r.t. total water matric potential - dPsiLiq_dPsi0 => deriv_data%var(iLookDERIV%dPsiLiq_dPsi0 )%dat ,& ! intent(out): [dp(:)] derivative in liquid water matric pot w.r.t. the total water matric pot (-) - dPsiLiq_dTemp => deriv_data%var(iLookDERIV%dPsiLiq_dTemp )%dat ,& ! intent(out): [dp(:)] derivative in the liquid water matric potential w.r.t. temperature - mLayerdTheta_dTk => deriv_data%var(iLookDERIV%mLayerdTheta_dTk)%dat ,& ! intent(out): [dp(:)] derivative of volumetric liquid water content w.r.t. temperature - dTheta_dTkCanopy => deriv_data%var(iLookDERIV%dTheta_dTkCanopy)%dat(1) ,& ! intent(out): [dp] derivative of volumetric liquid water content w.r.t. temperature - d2VolTot_d2Psi0 => deriv_data%var(iLookDERIV%d2VolTot_d2Psi0 )%dat ,& ! intent(out): [dp(:)] second derivative in total water content w.r.t. total water matric potential - mLayerd2Theta_dTk2 => deriv_data%var(iLookDERIV%mLayerd2Theta_dTk2 )%dat ,& ! intent(out): [dp(:)] second derivative of volumetric liquid water content w.r.t. temperature - d2Theta_dTkCanopy2 => deriv_data%var(iLookDERIV%d2Theta_dTkCanopy2 )%dat(1) ,& ! intent(out): [dp ] second derivative of volumetric liquid water content w.r.t. temperature - dFracLiqSnow_dTk => deriv_data%var(iLookDERIV%dFracLiqSnow_dTk )%dat ,& ! intent(out): [dp(:)] derivative in fraction of liquid snow w.r.t. temperature - dFracLiqVeg_dTkCanopy => deriv_data%var(iLookDERIV%dFracLiqVeg_dTkCanopy )%dat(1) ,& ! intent(out): [dp ] derivative in fraction of (throughfall + drainage) w.r.t. temperature - dVolHtCapBulk_dPsi0 => deriv_data%var(iLookDERIV%dVolHtCapBulk_dPsi0 )%dat ,& ! intent(out): [dp(:)] derivative in bulk heat capacity w.r.t. matric potential - dVolHtCapBulk_dTheta => deriv_data%var(iLookDERIV%dVolHtCapBulk_dTheta )%dat ,& ! intent(out): [dp(:)] derivative in bulk heat capacity w.r.t. volumetric water content - dVolHtCapBulk_dCanWat => deriv_data%var(iLookDERIV%dVolHtCapBulk_dCanWat)%dat(1) ,& ! intent(out): [dp ] derivative in bulk heat capacity w.r.t. volumetric water content - dVolHtCapBulk_dTk => deriv_data%var(iLookDERIV%dVolHtCapBulk_dTk )%dat ,& ! intent(out): [dp(:)] derivative in bulk heat capacity w.r.t. temperature - dVolHtCapBulk_dTkCanopy => deriv_data%var(iLookDERIV%dVolHtCapBulk_dTkCanopy)%dat(1) & ! intent(out): [dp ] derivative in bulk heat capacity w.r.t. temperature -) ! association with variables in the data structures - - ! -------------------------------------------------------------------------------------------------------------------------------- - ! -------------------------------------------------------------------------------------------------------------------------------- - - ! initialize error control - err=0; message='updateVars4JacDAE/' - - ! allocate space and assign values to the flag vector - allocate(computedCoupling(size(ixMapSubset2Full)),stat=err) ! .true. if computed the coupling for a given state variable - if(err/=0)then; message=trim(message)//'problem allocating computedCoupling'; return; endif - computedCoupling(:)=.false. - - ! loop through model state variables - do iState=1,size(ixMapSubset2Full) - - ! check the need for the computations - if(computedCoupling(iState)) cycle - - ! ----- - ! - compute indices... - ! -------------------- - - ! get domain type, and index of the control volume within the domain - ixFullVector = ixMapSubset2Full(iState) ! index within full state vector - ixDomainType = ixDomainType_subset(iState) ! named variables defining the domain (iname_cas, iname_veg, etc.) - ixControlIndex = ixControlVolume(ixFullVector) ! index within a given domain - - ! get the layer index - select case(ixDomainType) - case(iname_cas); cycle ! canopy air space: do nothing - case(iname_veg); iLayer = 0 - case(iname_snow); iLayer = ixControlIndex - case(iname_soil); iLayer = ixControlIndex + nSnow - case(iname_aquifer); cycle ! aquifer: do nothing - case default; err=20; message=trim(message)//'expect case to be iname_cas, iname_veg, iname_snow, iname_soil, iname_aquifer'; return - end select - - ! get the index of the other (energy or mass) state variable within the full state vector - select case(ixDomainType) - case(iname_veg) ; ixOther = merge(ixHydCanopy(1), ixNrgCanopy(1), ixStateType(ixFullVector)==iname_nrgCanopy) - case(iname_snow, iname_soil); ixOther = merge(ixHydLayer(iLayer),ixNrgLayer(iLayer),ixStateType(ixFullVector)==iname_nrgLayer) - case default; err=20; message=trim(message)//'expect case to be iname_veg, iname_snow, iname_soil'; return - end select - - ! get the index in the local state vector - ixOtherLocal = ixMapFull2Subset(ixOther) ! ixOtherLocal could equal integerMissing - if(ixOtherLocal/=integerMissing) computedCoupling(ixOtherLocal)=.true. - - ! check if we have a coupled solution - isCoupled = (ixOtherLocal/=integerMissing) - - ! check if we are an energy state - isNrgState = (ixStateType(ixFullVector)==iname_nrgCanopy .or. ixStateType(ixFullVector)==iname_nrgLayer) - - if(printFlag)then - print*, 'iState = ', iState, size(ixMapSubset2Full) - print*, 'ixFullVector = ', ixFullVector - print*, 'ixDomainType = ', ixDomainType - print*, 'ixControlIndex = ', ixControlIndex - print*, 'ixOther = ', ixOther - print*, 'ixOtherLocal = ', ixOtherLocal - print*, 'do_adjustTemp = ', do_adjustTemp - print*, 'isCoupled = ', isCoupled - print*, 'isNrgState = ', isNrgState - endif - - - - ! ======================================================================================================================================= - ! ======================================================================================================================================= - ! ======================================================================================================================================= - ! ======================================================================================================================================= - ! ======================================================================================================================================= - ! ======================================================================================================================================= - - ! update hydrology state variables for the uncoupled solution - if(.not.isNrgState .and. .not.isCoupled)then - - ! update the total water from volumetric liquid water - if(ixStateType(ixFullVector)==iname_liqCanopy .or. ixStateType(ixFullVector)==iname_liqLayer)then - select case(ixDomainType) - case(iname_veg) - scalarCanopyWatTrial = scalarCanopyLiqTrial + scalarCanopyIceTrial - scalarCanopyWatPrime = scalarCanopyLiqPrime + scalarCanopyIcePrime - case(iname_snow) - mLayerVolFracWatTrial(iLayer) = mLayerVolFracLiqTrial(iLayer) + mLayerVolFracIceTrial(iLayer)*iden_ice/iden_water - mLayerVolFracWatPrime(iLayer) = mLayerVolFracLiqPrime(iLayer) + mLayerVolFracIcePrime(iLayer)*iden_ice/iden_water - case(iname_soil) - mLayerVolFracWatTrial(iLayer) = mLayerVolFracLiqTrial(iLayer) + mLayerVolFracIceTrial(iLayer) ! no volume expansion - mLayerVolFracWatPrime(iLayer) = mLayerVolFracLiqPrime(iLayer) + mLayerVolFracIcePrime(iLayer) - case default; err=20; message=trim(message)//'expect case to be iname_veg, iname_snow, or iname_soil'; return - end select - endif - - ! update the total water and the total water matric potential - if(ixDomainType==iname_soil)then - select case( ixStateType(ixFullVector) ) - ! --> update the total water from the liquid water matric potential - case(iname_lmpLayer) - - effSat = volFracLiq(mLayerMatricHeadLiqTrial(ixControlIndex),vGn_alpha(ixControlIndex),0._rkind,1._rkind,vGn_n(ixControlIndex),vGn_m(ixControlIndex)) ! effective saturation - avPore = theta_sat(ixControlIndex) - mLayerVolFracIceTrial(iLayer) - theta_res(ixControlIndex) ! available pore space - mLayerVolFracLiqTrial(iLayer) = effSat*avPore + theta_res(ixControlIndex) - mLayerVolFracWatTrial(iLayer) = mLayerVolFracLiqTrial(iLayer) + mLayerVolFracIceTrial(iLayer) ! no volume expansion - mLayerVolFracWatPrime(iLayer) = mLayerVolFracLiqPrime(iLayer) + mLayerVolFracIcePrime(iLayer) - mLayerMatricHeadTrial(ixControlIndex) = matricHead(mLayerVolFracWatTrial(iLayer),vGn_alpha(ixControlIndex),theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex)) - mLayerMatricHeadPrime(ixControlIndex) = dPsi_dTheta(mLayerVolFracWatTrial(iLayer),vGn_alpha(ixControlIndex),theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex)) * mLayerVolFracWatPrime(iLayer) - !write(*,'(a,1x,i4,1x,3(f20.10,1x))') 'mLayerVolFracLiqTrial(iLayer) 1 = ', iLayer, mLayerVolFracLiqTrial(iLayer), mLayerVolFracIceTrial(iLayer), mLayerVolFracWatTrial(iLayer) - ! --> update the total water from the total water matric potential - case(iname_matLayer) - - mLayerVolFracWatTrial(iLayer) = volFracLiq(mLayerMatricHeadTrial(ixControlIndex),vGn_alpha(ixControlIndex),theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex)) - mLayerVolFracWatPrime(iLayer) = dTheta_dPsi(mLayerMatricHeadTrial(ixControlIndex),vGn_alpha(ixControlIndex),theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex)) *mLayerMatricHeadPrime(ixControlIndex) - ! --> update the total water matric potential (assume already have mLayerVolFracWatTrial given block above) - case(iname_liqLayer, iname_watLayer) - - mLayerMatricHeadTrial(ixControlIndex) = matricHead(mLayerVolFracWatTrial(iLayer),vGn_alpha(ixControlIndex),theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex)) - mLayerMatricHeadPrime(ixControlIndex) = dPsi_dTheta(mLayerVolFracWatTrial(iLayer),vGn_alpha(ixControlIndex),theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex)) * mLayerVolFracWatPrime(iLayer) - case default; err=20; message=trim(message)//'expect iname_lmpLayer, iname_matLayer, iname_liqLayer, or iname_watLayer'; return - end select - endif ! if in the soil domain - - endif ! if hydrology state variable or uncoupled solution - - ! compute the critical soil temperature below which ice exists - select case(ixDomainType) - case(iname_veg, iname_snow); Tcrit = Tfreeze; - case(iname_soil); Tcrit = crit_soilT( mLayerMatricHeadTrial(ixControlIndex) ); - case default; err=20; message=trim(message)//'expect case to be iname_veg, iname_snow, iname_soil'; return - end select - - ! initialize temperature - select case(ixDomainType) - case(iname_veg); xTemp = scalarCanopyTempTrial - case(iname_snow, iname_soil); xTemp = mLayerTempTrial(iLayer) - case default; err=20; message=trim(message)//'expect case to be iname_veg, iname_snow, iname_soil'; return - end select - - ! define brackets for the root - ! NOTE: start with an enormous range; updated quickly in the iterations - tempMin = xTemp - 10._rkind - tempMax = xTemp + 10._rkind - - ! get iterations (set to maximum iterations if adjusting the temperature) - niter = merge(maxiter, 1, do_adjustTemp) - - ! iterate - iterations: do iter=1,niter - - ! restrict temperature - if(xTemp <= tempMin .or. xTemp >= tempMax)then - xTemp = 0.5_rkind*(tempMin + tempMax) ! new value - bFlag = .true. - else - bFlag = .false. - endif - - ! ----- - ! - compute derivatives... - ! ------------------------ - - ! compute the derivative in bulk heat capacity w.r.t. total water content or water matric potential (m-1) - ! compute the derivative in total water content w.r.t. total water matric potential (m-1) - ! NOTE 1: valid for frozen and unfrozen conditions - ! NOTE 2: for case "iname_lmpLayer", dVolTot_dPsi0 = dVolLiq_dPsi, dVolHtCapBulk_dPsi0 may be wrong - select case(ixDomainType) - case(iname_veg) - fLiq = fracLiquid(xTemp,snowfrz_scale) - dVolHtCapBulk_dCanWat = ( -Cp_ice*( fLiq-1._rkind ) + Cp_water*fLiq )/canopyDepth !this is iden_water/(iden_water*canopyDepth) - case(iname_snow) - fLiq = fracLiquid(xTemp,snowfrz_scale) - dVolHtCapBulk_dTheta(iLayer) = iden_water * ( -Cp_ice*( fLiq-1._rkind ) + Cp_water*fLiq ) + iden_air * ( ( fLiq-1._rkind )*iden_water/iden_ice - fLiq ) * Cp_air - case(iname_soil) - select case( ixStateType(ixFullVector) ) - case(iname_lmpLayer) - dVolTot_dPsi0(ixControlIndex) = dTheta_dPsi(mLayerMatricHeadLiqTrial(ixControlIndex),vGn_alpha(ixControlIndex),0._rkind,1._rkind,vGn_n(ixControlIndex),vGn_m(ixControlIndex))*avPore - d2VolTot_d2Psi0(ixControlIndex) = d2Theta_dPsi2(mLayerMatricHeadLiqTrial(ixControlIndex),vGn_alpha(ixControlIndex),0._rkind,1._rkind,vGn_n(ixControlIndex),vGn_m(ixControlIndex))*avPore - case default - dVolTot_dPsi0(ixControlIndex) = dTheta_dPsi(mLayerMatricHeadTrial(ixControlIndex),vGn_alpha(ixControlIndex),theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex)) - d2VolTot_d2Psi0(ixControlIndex) = d2Theta_dPsi2(mLayerMatricHeadTrial(ixControlIndex),vGn_alpha(ixControlIndex),theta_res(ixControlIndex),theta_sat(ixControlIndex),& - vGn_n(ixControlIndex),vGn_m(ixControlIndex)) - end select - ! dVolHtCapBulk_dPsi0 uses the derivative in water retention curve above critical temp w.r.t.state variable dVolTot_dPsi0 - dVolHtCapBulk_dTheta(iLayer) = realMissing ! do not use - if(xTemp< Tcrit) dVolHtCapBulk_dPsi0(ixControlIndex) = (iden_ice * Cp_ice - iden_air * Cp_air) * dVolTot_dPsi0(ixControlIndex) - if(xTemp>=Tcrit) dVolHtCapBulk_dPsi0(ixControlIndex) = (iden_water * Cp_water - iden_air * Cp_air) * dVolTot_dPsi0(ixControlIndex) - end select - - ! compute the derivative in liquid water content w.r.t. temperature - ! --> partially frozen: dependence of liquid water on temperature - ! compute the derivative in bulk heat capacity w.r.t. temperature - if(xTemp<Tcrit)then - select case(ixDomainType) - case(iname_veg) - fLiq = fracLiquid(xTemp,snowfrz_scale) - dFracLiqVeg_dTkCanopy = dFracLiq_dTk(xTemp,snowfrz_scale) - dTheta_dTkCanopy = dFracLiqVeg_dTkCanopy * scalarCanopyWatTrial/(iden_water*canopyDepth) - d2Theta_dTkCanopy2 = 2._rkind * snowfrz_scale**2._rkind * ( (Tfreeze - xTemp) * 2._rkind * fLiq * dFracLiqVeg_dTkCanopy - fLiq**2._rkind ) * scalarCanopyWatTrial/(iden_water*canopyDepth) - dVolHtCapBulk_dTkCanopy = iden_water * (-Cp_ice + Cp_water) * dTheta_dTkCanopy !same as snow but there is no derivative in air - case(iname_snow) - fLiq = fracLiquid(xTemp,snowfrz_scale) - dFracLiqSnow_dTk(iLayer) = dFracLiq_dTk(xTemp,snowfrz_scale) - mLayerdTheta_dTk(iLayer) = dFracLiqSnow_dTk(iLayer) * mLayerVolFracWatTrial(iLayer) - mLayerd2Theta_dTk2(iLayer) = 2._rkind * snowfrz_scale**2._rkind * ( (Tfreeze - xTemp) * 2._rkind * fLiq * dFracLiqSnow_dTk(iLayer) - fLiq**2._rkind ) * mLayerVolFracWatTrial(iLayer) - dVolHtCapBulk_dTk(iLayer) = ( iden_water * (-Cp_ice + Cp_water) + iden_air * (iden_water/iden_ice - 1._rkind) * Cp_air ) * mLayerdTheta_dTk(iLayer) - case(iname_soil) - dFracLiqSnow_dTk(iLayer) = 0._rkind !dTheta_dTk(xTemp,theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_alpha(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex))/ mLayerVolFracWatTrial(iLayer) - mLayerdTheta_dTk(iLayer) = dTheta_dTk(xTemp,theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_alpha(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex)) - mLayerd2Theta_dTk2(iLayer) = d2Theta_dTk2(xTemp,theta_res(ixControlIndex),theta_sat(ixControlIndex),vGn_alpha(ixControlIndex),vGn_n(ixControlIndex),vGn_m(ixControlIndex)) - dVolHtCapBulk_dTk(iLayer) = (-iden_ice * Cp_ice + iden_water * Cp_water) * mLayerdTheta_dTk(iLayer) - case default; err=20; message=trim(message)//'expect case to be iname_veg, iname_snow, iname_soil'; return - end select ! domain type - - ! --> unfrozen: no dependence of liquid water on temperature - else - select case(ixDomainType) - case(iname_veg); dTheta_dTkCanopy = 0._rkind; d2Theta_dTkCanopy2 = 0._rkind; dFracLiqVeg_dTkCanopy = 0._rkind; dVolHtCapBulk_dTkCanopy = 0._rkind - case(iname_snow); mLayerdTheta_dTk(iLayer) = 0._rkind; mLayerd2Theta_dTk2(iLayer) = 0._rkind; dFracLiqSnow_dTk(iLayer) = 0._rkind; dVolHtCapBulk_dTk(iLayer) = 0._rkind - case(iname_soil); mLayerdTheta_dTk(iLayer) = 0._rkind; mLayerd2Theta_dTk2(iLayer) = 0._rkind; dFracLiqSnow_dTk(iLayer) = 0._rkind; dVolHtCapBulk_dTk(iLayer) = 0._rkind - end select ! domain type - endif - - ! ----- - ! - update volumetric fraction of liquid water and ice... - ! => case of hydrology state uncoupled with energy (and when not adjusting the temperature)... - ! ----------------------------------------------------------------------------------------------- - - ! case of hydrology state uncoupled with energy (and when not adjusting the temperature) - if(.not.do_adjustTemp .and. .not.isNrgState .and. .not.isCoupled)then - - ! compute the fraction of snow - select case(ixDomainType) - case(iname_veg); scalarFracLiqVeg = fracliquid(xTemp,snowfrz_scale) - case(iname_snow); mLayerFracLiqSnow(iLayer) = fracliquid(xTemp,snowfrz_scale) - case(iname_soil) ! do nothing - case default; err=20; message=trim(message)//'expect case to be iname_veg, iname_snow, iname_soil'; return - end select ! domain type - - ! ----- - ! - update volumetric fraction of liquid water and ice... - ! => case of energy state or coupled solution (or adjusting the temperature)... - ! -------------------------------------------------------------------------------- - - ! case of energy state OR coupled solution (or adjusting the temperature) - elseif(do_adjustTemp .or. ( (isNrgState .or. isCoupled) ) )then - - ! identify domain type - select case(ixDomainType) - - ! *** vegetation canopy - case(iname_veg) - ! compute volumetric fraction of liquid water and ice - call updateVegSundials(& - xTemp, & ! intent(in) : temperature (K) - scalarCanopyWatTrial, & ! intent(in) : canopy total water (kg m-2) - snowfrz_scale, & ! intent(in) : scaling parameter for the snow freezing curve (K-1) - scalarCanopyTempPrime, & ! intent(in) : canopy temperature time derivative (K/s) - scalarCanopyWatPrime, & ! intent(in) : canopy total water time derivative (kg m-2 /s) - scalarVolFracLiq, & ! intent(out) : trial canopy liquid water (-) - scalarVolFracIce, & ! intent(out) : trial volumetric canopy ice (-) - scalarVolFracLiqPrime, & ! intent(out) : trial volumetric canopy liquid water (-) - scalarVolFracIcePrime, & ! intent(out) : trial volumetric canopy ice (-) - scalarFracLiqVeg, & ! intent(out) : fraction of liquid water (-) - err,cmessage) ! intent(out) : error control - if(err/=0)then; message=trim(message)//trim(cmessage); return; endif - - ! *** snow layers - case(iname_snow) - - call updateSnowSundials(& - xTemp, & ! intent(in) : temperature (K) - mLayerVolFracWatTrial(iLayer), & ! intent(in) : mass state variable = trial volumetric fraction of water (-) - snowfrz_scale, & ! intent(in) : scaling parameter for the snow freezing curve (K-1) - mLayerTempPrime(iLayer), & ! intent(in) : temperature time derivative (K/s) - mLayerVolFracWatPrime(iLayer), & ! intent(in) : volumetric fraction of total water time derivative (-) - mLayerVolFracLiqTrial(iLayer), & ! intent(out) : trial volumetric fraction of liquid water (-) - mLayerVolFracIceTrial(iLayer), & ! intent(out) : trial volumetric fraction if ice (-) - mLayerVolFracLiqPrime(iLayer), & ! intent(out) : volumetric fraction of liquid water time derivative (-) - mLayerVolFracIcePrime(iLayer), & ! intent(out) : volumetric fraction of ice time derivative (-) - mLayerFracLiqSnow(iLayer), & ! intent(out) : fraction of liquid water (-) - err,cmessage) ! intent(out) : error control - if(err/=0)then; message=trim(message)//trim(cmessage); return; endif - - ! *** soil layers - case(iname_soil) - - ! compute volumetric fraction of liquid water and ice - call updateSoilSundials2(& - xTemp, & ! intent(in) : temperature (K) - mLayerMatricHeadTrial(ixControlIndex), & ! intent(in) : total water matric potential (m) - mLayerTempPrime(iLayer), & ! intent(in) : temperature time derivative (K/s) - mLayerMatricHeadPrime(ixControlIndex), & ! intent(in) : total water matric potential time derivative (m/s) - ! intent(in) : soil parameters - vGn_alpha(ixControlIndex), & ! intent(in) : van Genutchen "alpha" parameter - vGn_n(ixControlIndex), & ! intent(in) : van Genutchen "n" parameter - theta_sat(ixControlIndex), & ! intent(in) : soil porosity (-) - theta_res(ixControlIndex), & ! intent(in) : soil residual volumetric water content (-) - vGn_m(ixControlIndex), & ! intent(in) : van Genutchen "m" parameter (-) - mLayerVolFracWatTrial(iLayer), & ! intent(in) : mass state variable = trial volumetric fraction of water (-) - mLayerVolFracLiqTrial(iLayer), & ! intent(out) : trial volumetric fraction of liquid water (-) - mLayerVolFracIceTrial(iLayer), & ! intent(out) : trial volumetric fraction if ice (-) - mLayerVolFracWatPrime(iLayer), & ! intent(out) : volumetric fraction of total water time derivative (-) - mLayerVolFracLiqPrime(iLayer), & ! intent(out) : volumetric fraction of liquid water time derivative (-) - mLayerVolFracIcePrime(iLayer), & ! intent(out) : volumetric fraction of ice time derivative (-) - err,cmessage) ! intent(out) : error control - if(err/=0)then; message=trim(message)//trim(cmessage); return; endif - - ! check - case default; err=20; message=trim(message)//'expect case to be iname_veg, iname_snow, iname_soil'; return - - end select ! domain type - - ! final check - else - - ! do nothing (input = output) -- and check that we got here correctly - if( (isNrgState .or. isCoupled) )then - scalarVolFracLiq = realMissing - scalarVolFracIce = realMissing - else - message=trim(message)//'unexpected else branch' - err=20; return - endif - - endif ! if energy state or solution is coupled - - ! ----- - ! - update temperatures... - ! ------------------------ - - ! check the need to adjust temperature - if(do_adjustTemp)then - - ! get the melt energy - meltNrg = merge(LH_fus*iden_ice, LH_fus*iden_water, ixDomainType==iname_snow) - - ! compute the residual and the derivative - select case(ixDomainType) - - ! * vegetation - case(iname_veg) - call xTempSolve(& - ! constant over iterations - meltNrg = meltNrg ,& ! intent(in) : energy for melt+freeze (J m-3) - heatCap = scalarBulkVolHeatCapVeg ,& ! intent(in) : volumetric heat capacity (J m-3 K-1) - tempInit = scalarCanopyTemp ,& ! intent(in) : initial temperature (K) - volFracIceInit = scalarCanopyIce/(iden_water*canopyDepth),& ! intent(in) : initial volumetric fraction of ice (-) - ! trial values - xTemp = xTemp ,& ! intent(inout) : trial value of temperature - dLiq_dT = dTheta_dTkCanopy ,& ! intent(in) : derivative in liquid water content w.r.t. temperature (K-1) - volFracIceTrial = scalarVolFracIce ,& ! intent(in) : trial value for volumetric fraction of ice - ! residual and derivative - residual = residual ,& ! intent(out) : residual (J m-3) - derivative = derivative ) ! intent(out) : derivative (J m-3 K-1) - - ! * snow and soil - case(iname_snow, iname_soil) - call xTempSolve(& - ! constant over iterations - meltNrg = meltNrg ,& ! intent(in) : energy for melt+freeze (J m-3) - heatCap = mLayerVolHtCapBulk(iLayer) ,& ! intent(in) : volumetric heat capacity (J m-3 K-1) - tempInit = mLayerTemp(iLayer) ,& ! intent(in) : initial temperature (K) - volFracIceInit = mLayerVolFracIce(iLayer) ,& ! intent(in) : initial volumetric fraction of ice (-) - ! trial values - xTemp = xTemp ,& ! intent(inout) : trial value of temperature - dLiq_dT = mLayerdTheta_dTk(iLayer) ,& ! intent(in) : derivative in liquid water content w.r.t. temperature (K-1) - volFracIceTrial = mLayerVolFracIceTrial(iLayer) ,& ! intent(in) : trial value for volumetric fraction of ice - ! residual and derivative - residual = residual ,& ! intent(out) : residual (J m-3) - derivative = derivative ) ! intent(out) : derivative (J m-3 K-1) - - ! * check - case default; err=20; message=trim(message)//'expect case to be iname_veg, iname_snow, iname_soil'; return - - end select ! domain type - - ! check validity of residual - if( ieee_is_nan(residual) )then - message=trim(message)//'residual is not valid' - err=20; return - endif - - ! update bracket - if(residual < 0._rkind)then - tempMax = min(xTemp,tempMax) - else - tempMin = max(tempMin,xTemp) - end if - - ! compute iteration increment - tempInc = residual/derivative ! K - - ! check - if(globalPrintFlag)& - write(*,'(i4,1x,e20.10,1x,5(f20.10,1x),L1)') iter, residual, xTemp-Tcrit, tempInc, Tcrit, tempMin, tempMax, bFlag - - ! check convergence - if(abs(residual) < nrgConvTol .or. abs(tempInc) < tempConvTol) exit iterations - - ! add constraints for snow temperature - if(ixDomainType==iname_veg .or. ixDomainType==iname_snow)then - if(tempInc > Tcrit - xTemp) tempInc=(Tcrit - xTemp)*0.5_rkind ! simple bi-section method - endif ! if the domain is vegetation or snow - - ! deal with the discontinuity between partially frozen and unfrozen soil - if(ixDomainType==iname_soil)then - ! difference from the temperature below which ice exists - critDiff = Tcrit - xTemp - ! --> initially frozen (T < Tcrit) - if(critDiff > 0._rkind)then - if(tempInc > critDiff) tempInc = critDiff + epsT ! set iteration increment to slightly above critical temperature - ! --> initially unfrozen (T > Tcrit) - else - if(tempInc < critDiff) tempInc = critDiff - epsT ! set iteration increment to slightly below critical temperature - endif - endif ! if the domain is soil - - ! update the temperature trial - xTemp = xTemp + tempInc - - ! check failed convergence - if(iter==maxiter)then - message=trim(message)//'failed to converge' - err=-20; return ! negative error code = try to recover - endif - - endif ! if adjusting the temperature - - end do iterations ! iterating - - ! save temperature - select case(ixDomainType) - case(iname_veg); scalarCanopyTempTrial = xTemp - case(iname_snow, iname_soil); mLayerTempTrial(iLayer) = xTemp - end select - - ! ======================================================================================================================================= - ! ======================================================================================================================================= - - ! ----- - ! - compute the liquid water matric potential (and necessay derivatives)... - ! ------------------------------------------------------------------------- - - ! only for soil - if(ixDomainType==iname_soil)then - - ! check liquid water - if(mLayerVolFracLiqTrial(iLayer) > theta_sat(ixControlIndex) )then - message=trim(message)//'liquid water greater than porosity' - err=20; return - endif - - ! case of hydrology state uncoupled with energy - if(.not.isNrgState .and. .not.isCoupled)then - - ! derivatives relating liquid water matric potential to total water matric potential and temperature - dPsiLiq_dPsi0(ixControlIndex) = 1._rkind ! exact correspondence (psiLiq=psi0) - dPsiLiq_dTemp(ixControlIndex) = 0._rkind ! no relationship between liquid water matric potential and temperature - - ! case of energy state or coupled solution - else - ! compute the liquid matric potential (and the derivatives w.r.t. total matric potential and temperature) - call liquidHeadSundials(& - ! input - mLayerMatricHeadTrial(ixControlIndex) ,& ! intent(in) : total water matric potential (m) - mLayerMatricHeadPrime(ixControlIndex) ,& ! - mLayerVolFracLiqTrial(iLayer) ,& ! intent(in) : volumetric fraction of liquid water (-) - mLayerVolFracIceTrial(iLayer) ,& ! intent(in) : volumetric fraction of ice (-) - vGn_alpha(ixControlIndex),vGn_n(ixControlIndex),theta_sat(ixControlIndex),theta_res(ixControlIndex),vGn_m(ixControlIndex), & ! intent(in) : soil parameters - dVolTot_dPsi0(ixControlIndex) ,& ! intent(in) : derivative in the soil water characteristic (m-1) - mLayerdTheta_dTk(iLayer) ,& ! intent(in) : derivative in volumetric total water w.r.t. temperature (K-1) - mLayerVolFracLiqPrime(iLayer) ,& - mLayerVolFracIcePrime(iLayer) ,& - ! output - mLayerMatricHeadLiqTrial(ixControlIndex) ,& ! intent(out): liquid water matric potential (m) - mLayerMatricHeadLiqPrime(ixControlIndex) ,& ! - dPsiLiq_dPsi0(ixControlIndex) ,& ! intent(out): derivative in the liquid water matric potential w.r.t. the total water matric potential (-) - dPsiLiq_dTemp(ixControlIndex) ,& ! intent(out): derivative in the liquid water matric potential w.r.t. temperature (m K-1) - err,cmessage) ! intent(out): error control - if(err/=0)then; message=trim(message)//trim(cmessage); return; endif - - endif ! switch between hydrology and energy state - - endif ! if domain is soil - - end do ! looping through state variables - - ! deallocate space - deallocate(computedCoupling,stat=err) ! .true. if computed the coupling for a given state variable - if(err/=0)then; message=trim(message)//'problem deallocating computedCoupling'; return; endif - - ! end association to the variables in the data structures - end associate - - end subroutine updateVars4JacDAE - - - ! ********************************************************************************************************** - ! private subroutine xTempSolve: compute residual and derivative for temperature - ! ********************************************************************************************************** - subroutine xTempSolve(& - ! input: constant over iterations - meltNrg ,& ! intent(in) : energy for melt+freeze (J m-3) - heatCap ,& ! intent(in) : volumetric heat capacity (J m-3 K-1) - tempInit ,& ! intent(in) : initial temperature (K) - volFracIceInit ,& ! intent(in) : initial volumetric fraction of ice (-) - ! input-output: trial values - xTemp ,& ! intent(inout) : trial value of temperature - dLiq_dT ,& ! intent(in) : derivative in liquid water content w.r.t. temperature (K-1) - volFracIceTrial ,& ! intent(in) : trial value for volumetric fraction of ice - ! output: residual and derivative - residual ,& ! intent(out) : residual (J m-3) - derivative ) ! intent(out) : derivative (J m-3 K-1) - implicit none - ! input: constant over iterations - real(rkind),intent(in) :: meltNrg ! energy for melt+freeze (J m-3) - real(rkind),intent(in) :: heatCap ! volumetric heat capacity (J m-3 K-1) - real(rkind),intent(in) :: tempInit ! initial temperature (K) - real(rkind),intent(in) :: volFracIceInit ! initial volumetric fraction of ice (-) - ! input-output: trial values - real(rkind),intent(inout) :: xTemp ! trial value for temperature - real(rkind),intent(in) :: dLiq_dT ! derivative in liquid water content w.r.t. temperature (K-1) - real(rkind),intent(in) :: volFracIceTrial ! trial value for the volumetric fraction of ice (-) - ! output: residual and derivative - real(rkind),intent(out) :: residual ! residual (J m-3) - real(rkind),intent(out) :: derivative ! derivative (J m-3 K-1) - ! subroutine starts here - residual = -heatCap*(xTemp - tempInit) + meltNrg*(volFracIceTrial - volFracIceInit) ! J m-3 - derivative = heatCap + LH_fus*iden_water*dLiq_dT ! J m-3 K-1 - end subroutine xTempSolve - -end module updateVars4JacDAE_module diff --git a/utils/laugh_tests/celia1990/verification_data/runinfo.txt b/utils/laugh_tests/celia1990/verification_data/runinfo.txt index 5a1f0568f8e559fff7e56f1fc7d2cbe84308981d..66a332f2d5c81a4b4894f0651e81cf2a28d0d17d 100644 --- a/utils/laugh_tests/celia1990/verification_data/runinfo.txt +++ b/utils/laugh_tests/celia1990/verification_data/runinfo.txt @@ -1 +1 @@ - Run start time on system: ccyy=2022 - mm=09 - dd=09 - hh=14 - mi=49 - ss=55.809 + Run start time on system: ccyy=2022 - mm=09 - dd=09 - hh=20 - mi=11 - ss=42.068 diff --git a/utils/laugh_tests/celia1990/verification_data/summa_celia1990_G1-1_timestep.nc b/utils/laugh_tests/celia1990/verification_data/summa_celia1990_G1-1_timestep.nc index 87a123f7d4f6fa9f9397871cc5b572de51cb794b..3d0fbbe5494ae0c1443aedfd3f198532f0193d8b 100644 Binary files a/utils/laugh_tests/celia1990/verification_data/summa_celia1990_G1-1_timestep.nc and b/utils/laugh_tests/celia1990/verification_data/summa_celia1990_G1-1_timestep.nc differ diff --git a/utils/laugh_tests/colbeck1976/verification_data/runinfo.txt b/utils/laugh_tests/colbeck1976/verification_data/runinfo.txt index 537c18b416ed451203214dc90a55dd2a2b604d10..25874f090a496bf205b34945ebc0f70159fdebef 100644 --- a/utils/laugh_tests/colbeck1976/verification_data/runinfo.txt +++ b/utils/laugh_tests/colbeck1976/verification_data/runinfo.txt @@ -1 +1 @@ - Run start time on system: ccyy=2022 - mm=09 - dd=09 - hh=18 - mi=08 - ss=14.611 + Run start time on system: ccyy=2022 - mm=09 - dd=09 - hh=20 - mi=06 - ss=06.185 diff --git a/utils/laugh_tests/mizoguchi1990/verification_data/mizoguchi1990_G1-1_timestep.nc b/utils/laugh_tests/mizoguchi1990/verification_data/mizoguchi1990_G1-1_timestep.nc index 8afdcaad6cd4d3e56cb818cea27fcdbdec51a7ea..1df5a2177fde0cf94afe8c1146adc43674caac48 100644 Binary files a/utils/laugh_tests/mizoguchi1990/verification_data/mizoguchi1990_G1-1_timestep.nc and b/utils/laugh_tests/mizoguchi1990/verification_data/mizoguchi1990_G1-1_timestep.nc differ diff --git a/utils/laugh_tests/mizoguchi1990/verification_data/runinfo.txt b/utils/laugh_tests/mizoguchi1990/verification_data/runinfo.txt index 9d14b1d1e2194bd1de9dfd130c365ed5083147ae..63c3ff58e0f2b2fcbdca1b1c9de5143160b3aa13 100644 --- a/utils/laugh_tests/mizoguchi1990/verification_data/runinfo.txt +++ b/utils/laugh_tests/mizoguchi1990/verification_data/runinfo.txt @@ -1 +1 @@ - Run start time on system: ccyy=2022 - mm=09 - dd=09 - hh=16 - mi=18 - ss=28.178 + Run start time on system: ccyy=2022 - mm=09 - dd=09 - hh=20 - mi=14 - ss=07.920