From 054a97c92d105229114265a5fb8b4ec81760f2d2 Mon Sep 17 00:00:00 2001
From: Kyle <kyle.c.klenk@gmail.com>
Date: Mon, 29 Aug 2022 18:46:09 +0000
Subject: [PATCH] added error messeges to coupled em
---
build/source/driver/SummaActors_modelRun.f90 | 17 +-
build/source/engine/coupled_em.f90 | 1484 ++++++++++--------
build/source/engine/nr_utility.f90 | 17 -
3 files changed, 861 insertions(+), 657 deletions(-)
diff --git a/build/source/driver/SummaActors_modelRun.f90 b/build/source/driver/SummaActors_modelRun.f90
index 3310197..d5a48b3 100755
--- a/build/source/driver/SummaActors_modelRun.f90
+++ b/build/source/driver/SummaActors_modelRun.f90
@@ -187,7 +187,11 @@ contains
fracJulDay, &
yearLength, &
err,cmessage) ! intent(out): error control
- if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
+ if(err/=0)then
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ endif
! save the flag for computing the vegetation fluxes
if(computeVegFluxFlag) computeVegFlux = yes
@@ -215,7 +219,7 @@ contains
bvarStruct%var(iLookBVAR%basin__SurfaceRunoff)%dat(1) = 0._dp ! surface runoff (m s-1)
bvarStruct%var(iLookBVAR%basin__SoilDrainage)%dat(1) = 0._dp
bvarStruct%var(iLookBVAR%basin__ColumnOutflow)%dat(1) = 0._dp ! outflow from all "outlet" HRUs (those with no downstream HRU)
- bvarStruct%var(iLookBVAR%basin__TotalRunoff)%dat(1) = 0._dp
+ bvarStruct%var(iLookBVAR%basin__TotalRunoff)%dat(1) = 0._dp
! initialize baseflow variables
bvarStruct%var(iLookBVAR%basin__AquiferRecharge)%dat(1) = 0._dp ! recharge to the aquifer (m s-1)
@@ -244,6 +248,7 @@ contains
allocate(zSoilReverseSign(nSoil),stat=err)
if(err/=0)then
message=trim(message)//'problem allocating space for zSoilReverseSign'
+ print*, message
err=20; return
endif
zSoilReverseSign(:) = -progStruct%var(iLookPROG%iLayerHeight)%dat(nSnow+1:nLayers)
@@ -262,6 +267,7 @@ contains
deallocate(zSoilReverseSign,stat=err)
if(err/=0)then
message=trim(message)//'problem deallocating space for zSoilReverseSign'
+ print*, message
err=20; return
endif
@@ -290,7 +296,12 @@ contains
fluxStruct, & ! data structure of model fluxes
tmZoneOffsetFracDay,&
err,cmessage) ! error control
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
+ if(err/=0)then
+ err=20
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+endif
! initialize the number of flux calls
diagStruct%var(iLookDIAG%numFluxCalls)%dat(1) = 0._dp
diff --git a/build/source/engine/coupled_em.f90 b/build/source/engine/coupled_em.f90
index de184bf..964bb2e 100755
--- a/build/source/engine/coupled_em.f90
+++ b/build/source/engine/coupled_em.f90
@@ -120,333 +120,371 @@ subroutine coupled_em(&
yearLength, &
! error control
err,message) ! intent(out): error control
- ! structure allocations
- USE allocspace4chm_module,only:allocLocal ! allocate local data structures
- USE allocspace4chm_module,only:resizeData ! clone a data structure
- ! preliminary subroutines
- USE vegPhenlgy_module,only:vegPhenlgy ! compute vegetation phenology
- USE vegNrgFlux_module,only:wettedFrac ! compute wetted fraction of the canopy (used in sw radiation fluxes)
- USE snowAlbedo_module,only:snowAlbedo ! compute snow albedo
- USE vegSWavRad_module,only:vegSWavRad ! compute canopy sw radiation fluxes
- USE canopySnow_module,only:canopySnow ! compute interception and unloading of snow from the vegetation canopy
- USE volicePack_module,only:newsnwfall ! compute change in the top snow layer due to throughfall and unloading
- USE volicePack_module,only:volicePack ! merge and sub-divide snow layers, if necessary
- USE diagn_evar_module,only:diagn_evar ! compute diagnostic energy variables -- thermal conductivity and heat capacity
- ! the model solver
- USE indexState_module,only:indexState ! define indices for all model state variables and layers
- USE opSplittin_module,only:opSplittin ! solve the system of thermodynamic and hydrology equations for a given substep
- ! additional subroutines
- USE tempAdjust_module,only:tempAdjust ! adjust snow temperature associated with new snowfall
- USE snwDensify_module,only:snwDensify ! snow densification (compaction and cavitation)
- USE var_derive_module,only:calcHeight ! module to calculate height at layer interfaces and layer mid-point
- ! look-up values for the numerical method
- USE mDecisions_module,only: &
- iterative, & ! iterative
- nonIterative, & ! non-iterative
- iterSurfEnergyBal ! iterate only on the surface energy balance
- ! look-up values for the maximum interception capacity
- USE mDecisions_module,only: &
- stickySnow, & ! maximum interception capacity an increasing function of temerature
- lightSnow ! maximum interception capacity an inverse function of new snow density
- implicit none
- ! model control
- integer(4),intent(in) :: indxHRU ! hruId
- real(dp),intent(inout) :: dt_init ! used to initialize the size of the sub-step
- integer(i4b),intent(in) :: dt_init_factor ! Used to adjust the length of the timestep in the event of a failure
- logical(lgt),intent(inout) :: computeVegFlux ! flag to indicate if we are computing fluxes over vegetation (.false. means veg is buried with snow)
- ! data structures (input)
- type(var_i),intent(in) :: type_data ! type of vegetation and soil
- type(var_d),intent(in) :: attr_data ! spatial attributes
- type(var_d),intent(in) :: forc_data ! model forcing data
- type(var_dlength),intent(in) :: mpar_data ! model parameters
- type(var_dlength),intent(in) :: bvar_data ! basin-average model variables
- ! data structures (input-output)
- type(var_ilength),intent(inout) :: indx_data ! state vector geometry
- type(var_dlength),intent(inout) :: 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) :: flux_data ! model fluxes for a local HRU
- real(dp),intent(inout) :: fracJulDay
- integer(i4b),intent(inout) :: yearLength
- ! error control
- integer(i4b),intent(out) :: err ! error code
- character(*),intent(out) :: message ! error message
- ! =====================================================================================================================================================
- ! =====================================================================================================================================================
- ! local variables
- character(len=256) :: cmessage ! error message
- integer(i4b) :: nSnow ! number of snow layers
- integer(i4b) :: nSoil ! number of soil layers
- integer(i4b) :: nLayers ! total number of layers
- integer(i4b) :: nState ! total number of state variables
- real(dp) :: dtSave ! length of last input model sub-step (seconds)
- real(dp) :: dt_sub ! length of model sub-step (seconds)
- real(dp) :: dt_wght ! weight applied to model sub-step (dt_sub/data_step)
- real(dp) :: dt_solv ! seconds in the data step that have been completed
- real(dp) :: dtMultiplier ! time step multiplier (-) based on what happenned in "opSplittin"
- real(dp) :: minstep,maxstep ! minimum and maximum time step length (seconds)
- integer(i4b) :: nsub ! number of substeps
- logical(lgt) :: computeVegFluxOld ! flag to indicate if we are computing fluxes over vegetation on the previous sub step
- logical(lgt) :: includeAquifer ! flag to denote that an aquifer is included
- logical(lgt) :: modifiedLayers ! flag to denote that snow layers were modified
- logical(lgt) :: modifiedVegState ! flag to denote that vegetation states were modified
- type(var_dlength) :: flux_mean ! timestep-average model fluxes for a local HRU
- integer(i4b) :: nLayersRoots ! number of soil layers that contain roots
- real(dp) :: exposedVAI ! exposed vegetation area index
- real(dp) :: dCanopyWetFraction_dWat ! derivative in wetted fraction w.r.t. canopy total water (kg-1 m2)
- real(dp) :: dCanopyWetFraction_dT ! derivative in wetted fraction w.r.t. canopy temperature (K-1)
- real(dp),parameter :: varNotUsed1=-9999._dp ! variables used to calculate derivatives (not needed here)
- real(dp),parameter :: varNotUsed2=-9999._dp ! variables used to calculate derivatives (not needed here)
- integer(i4b) :: iSnow ! index of snow layers
- integer(i4b) :: iLayer ! index of model layers
- real(dp) :: massLiquid ! mass liquid water (kg m-2)
- real(dp) :: superflousSub ! superflous sublimation (kg m-2 s-1)
- real(dp) :: superflousNrg ! superflous energy that cannot be used for sublimation (W m-2 [J m-2 s-1])
- integer(i4b) :: ixSolution ! solution method used by opSplitting
- logical(lgt) :: firstSubStep ! flag to denote if the first time step
- logical(lgt) :: stepFailure ! flag to denote the need to reduce length of the coupled step and try again
- logical(lgt) :: tooMuchMelt ! flag to denote that there was too much melt in a given time step
- logical(lgt) :: doLayerMerge ! flag to denote the need to merge snow layers
- logical(lgt) :: pauseFlag ! flag to pause execution
- logical(lgt),parameter :: backwardsCompatibility=.true. ! flag to denote a desire to ensure backwards compatibility with previous branches.
- type(var_ilength) :: indx_temp ! temporary model index variables
- type(var_dlength) :: prog_temp ! temporary model prognostic variables
- type(var_dlength) :: diag_temp ! temporary model diagnostic variables
- ! check SWE
- real(dp) :: oldSWE ! SWE at the start of the substep
- real(dp) :: newSWE ! SWE at the end of the substep
- real(dp) :: delSWE ! change in SWE over the subtep
- real(dp) :: effRainfall ! effective rainfall (kg m-2 s-1)
- real(dp) :: effSnowfall ! effective snowfall (kg m-2 s-1)
- real(dp) :: sfcMeltPond ! surface melt pond (kg m-2)
- real(dp) :: massBalance ! mass balance error (kg m-2)
- ! balance checks
- integer(i4b) :: iVar ! loop through model variables
- real(dp) :: totalSoilCompress ! total soil compression (kg m-2)
- real(dp) :: scalarCanopyWatBalError ! water balance error for the vegetation canopy (kg m-2)
- real(dp) :: scalarSoilWatBalError ! water balance error (kg m-2)
- real(dp) :: scalarInitCanopyLiq ! initial liquid water on the vegetation canopy (kg m-2)
- real(dp) :: scalarInitCanopyIce ! initial ice on the vegetation canopy (kg m-2)
- real(dp) :: balanceCanopyWater0 ! total water stored in the vegetation canopy at the start of the step (kg m-2)
- real(dp) :: balanceCanopyWater1 ! total water stored in the vegetation canopy at the end of the step (kg m-2)
- real(dp) :: balanceSoilWater0 ! total soil storage at the start of the step (kg m-2)
- real(dp) :: balanceSoilWater1 ! total soil storage at the end of the step (kg m-2)
- real(dp) :: balanceSoilInflux ! input to the soil zone
- real(dp) :: balanceSoilBaseflow ! output from the soil zone
- real(dp) :: balanceSoilDrainage ! output from the soil zone
- real(dp) :: balanceSoilET ! output from the soil zone
- real(dp) :: balanceAquifer0 ! total aquifer storage at the start of the step (kg m-2)
- real(dp) :: balanceAquifer1 ! total aquifer storage at the end of the step (kg m-2)
- ! test balance checks
- logical(lgt), parameter :: printBalance=.false. ! flag to print the balance checks
- real(dp), allocatable :: liqSnowInit(:) ! volumetric liquid water conetnt of snow at the start of the time step
- real(dp), allocatable :: liqSoilInit(:) ! soil moisture at the start of the time step
- ! ----------------------------------------------------------------------------------------------------------------------------------------------
- ! initialize error control
- err=0; message="coupled_em/"
+ ! structure allocations
+ USE allocspace4chm_module,only:allocLocal ! allocate local data structures
+ USE allocspace4chm_module,only:resizeData ! clone a data structure
+ ! preliminary subroutines
+ USE vegPhenlgy_module,only:vegPhenlgy ! compute vegetation phenology
+ USE vegNrgFlux_module,only:wettedFrac ! compute wetted fraction of the canopy (used in sw radiation fluxes)
+ USE snowAlbedo_module,only:snowAlbedo ! compute snow albedo
+ USE vegSWavRad_module,only:vegSWavRad ! compute canopy sw radiation fluxes
+ USE canopySnow_module,only:canopySnow ! compute interception and unloading of snow from the vegetation canopy
+ USE volicePack_module,only:newsnwfall ! compute change in the top snow layer due to throughfall and unloading
+ USE volicePack_module,only:volicePack ! merge and sub-divide snow layers, if necessary
+ USE diagn_evar_module,only:diagn_evar ! compute diagnostic energy variables -- thermal conductivity and heat capacity
+ ! the model solver
+ USE indexState_module,only:indexState ! define indices for all model state variables and layers
+ USE opSplittin_module,only:opSplittin ! solve the system of thermodynamic and hydrology equations for a given substep
+ ! additional subroutines
+ USE tempAdjust_module,only:tempAdjust ! adjust snow temperature associated with new snowfall
+ USE snwDensify_module,only:snwDensify ! snow densification (compaction and cavitation)
+ USE var_derive_module,only:calcHeight ! module to calculate height at layer interfaces and layer mid-point
+ ! look-up values for the numerical method
+ USE mDecisions_module,only: &
+ iterative, & ! iterative
+ nonIterative, & ! non-iterative
+ iterSurfEnergyBal ! iterate only on the surface energy balance
+ ! look-up values for the maximum interception capacity
+ USE mDecisions_module,only: &
+ stickySnow, & ! maximum interception capacity an increasing function of temerature
+ lightSnow ! maximum interception capacity an inverse function of new snow density
+ implicit none
+ ! model control
+ integer(4),intent(in) :: indxHRU ! hruId
+ real(dp),intent(inout) :: dt_init ! used to initialize the size of the sub-step
+ integer(i4b),intent(in) :: dt_init_factor ! Used to adjust the length of the timestep in the event of a failure
+ logical(lgt),intent(inout) :: computeVegFlux ! flag to indicate if we are computing fluxes over vegetation (.false. means veg is buried with snow)
+ ! data structures (input)
+ type(var_i),intent(in) :: type_data ! type of vegetation and soil
+ type(var_d),intent(in) :: attr_data ! spatial attributes
+ type(var_d),intent(in) :: forc_data ! model forcing data
+ type(var_dlength),intent(in) :: mpar_data ! model parameters
+ type(var_dlength),intent(in) :: bvar_data ! basin-average model variables
+ ! data structures (input-output)
+ type(var_ilength),intent(inout) :: indx_data ! state vector geometry
+ type(var_dlength),intent(inout) :: 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) :: flux_data ! model fluxes for a local HRU
+ real(dp),intent(inout) :: fracJulDay
+ integer(i4b),intent(inout) :: yearLength
+ ! error control
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! =====================================================================================================================================================
+ ! =====================================================================================================================================================
+ ! local variables
+ character(len=256) :: cmessage ! error message
+ integer(i4b) :: nSnow ! number of snow layers
+ integer(i4b) :: nSoil ! number of soil layers
+ integer(i4b) :: nLayers ! total number of layers
+ integer(i4b) :: nState ! total number of state variables
+ real(dp) :: dtSave ! length of last input model sub-step (seconds)
+ real(dp) :: dt_sub ! length of model sub-step (seconds)
+ real(dp) :: dt_wght ! weight applied to model sub-step (dt_sub/data_step)
+ real(dp) :: dt_solv ! seconds in the data step that have been completed
+ real(dp) :: dtMultiplier ! time step multiplier (-) based on what happenned in "opSplittin"
+ real(dp) :: minstep,maxstep ! minimum and maximum time step length (seconds)
+ integer(i4b) :: nsub ! number of substeps
+ logical(lgt) :: computeVegFluxOld ! flag to indicate if we are computing fluxes over vegetation on the previous sub step
+ logical(lgt) :: includeAquifer ! flag to denote that an aquifer is included
+ logical(lgt) :: modifiedLayers ! flag to denote that snow layers were modified
+ logical(lgt) :: modifiedVegState ! flag to denote that vegetation states were modified
+ type(var_dlength) :: flux_mean ! timestep-average model fluxes for a local HRU
+ integer(i4b) :: nLayersRoots ! number of soil layers that contain roots
+ real(dp) :: exposedVAI ! exposed vegetation area index
+ real(dp) :: dCanopyWetFraction_dWat ! derivative in wetted fraction w.r.t. canopy total water (kg-1 m2)
+ real(dp) :: dCanopyWetFraction_dT ! derivative in wetted fraction w.r.t. canopy temperature (K-1)
+ real(dp),parameter :: varNotUsed1=-9999._dp ! variables used to calculate derivatives (not needed here)
+ real(dp),parameter :: varNotUsed2=-9999._dp ! variables used to calculate derivatives (not needed here)
+ integer(i4b) :: iSnow ! index of snow layers
+ integer(i4b) :: iLayer ! index of model layers
+ real(dp) :: massLiquid ! mass liquid water (kg m-2)
+ real(dp) :: superflousSub ! superflous sublimation (kg m-2 s-1)
+ real(dp) :: superflousNrg ! superflous energy that cannot be used for sublimation (W m-2 [J m-2 s-1])
+ integer(i4b) :: ixSolution ! solution method used by opSplitting
+ logical(lgt) :: firstSubStep ! flag to denote if the first time step
+ logical(lgt) :: stepFailure ! flag to denote the need to reduce length of the coupled step and try again
+ logical(lgt) :: tooMuchMelt ! flag to denote that there was too much melt in a given time step
+ logical(lgt) :: doLayerMerge ! flag to denote the need to merge snow layers
+ logical(lgt) :: pauseFlag ! flag to pause execution
+ logical(lgt),parameter :: backwardsCompatibility=.true. ! flag to denote a desire to ensure backwards compatibility with previous branches.
+ type(var_ilength) :: indx_temp ! temporary model index variables
+ type(var_dlength) :: prog_temp ! temporary model prognostic variables
+ type(var_dlength) :: diag_temp ! temporary model diagnostic variables
+ ! check SWE
+ real(dp) :: oldSWE ! SWE at the start of the substep
+ real(dp) :: newSWE ! SWE at the end of the substep
+ real(dp) :: delSWE ! change in SWE over the subtep
+ real(dp) :: effRainfall ! effective rainfall (kg m-2 s-1)
+ real(dp) :: effSnowfall ! effective snowfall (kg m-2 s-1)
+ real(dp) :: sfcMeltPond ! surface melt pond (kg m-2)
+ real(dp) :: massBalance ! mass balance error (kg m-2)
+ ! balance checks
+ integer(i4b) :: iVar ! loop through model variables
+ real(dp) :: totalSoilCompress ! total soil compression (kg m-2)
+ real(dp) :: scalarCanopyWatBalError ! water balance error for the vegetation canopy (kg m-2)
+ real(dp) :: scalarSoilWatBalError ! water balance error (kg m-2)
+ real(dp) :: scalarInitCanopyLiq ! initial liquid water on the vegetation canopy (kg m-2)
+ real(dp) :: scalarInitCanopyIce ! initial ice on the vegetation canopy (kg m-2)
+ real(dp) :: balanceCanopyWater0 ! total water stored in the vegetation canopy at the start of the step (kg m-2)
+ real(dp) :: balanceCanopyWater1 ! total water stored in the vegetation canopy at the end of the step (kg m-2)
+ real(dp) :: balanceSoilWater0 ! total soil storage at the start of the step (kg m-2)
+ real(dp) :: balanceSoilWater1 ! total soil storage at the end of the step (kg m-2)
+ real(dp) :: balanceSoilInflux ! input to the soil zone
+ real(dp) :: balanceSoilBaseflow ! output from the soil zone
+ real(dp) :: balanceSoilDrainage ! output from the soil zone
+ real(dp) :: balanceSoilET ! output from the soil zone
+ real(dp) :: balanceAquifer0 ! total aquifer storage at the start of the step (kg m-2)
+ real(dp) :: balanceAquifer1 ! total aquifer storage at the end of the step (kg m-2)
+ ! test balance checks
+ logical(lgt), parameter :: printBalance=.false. ! flag to print the balance checks
+ real(dp), allocatable :: liqSnowInit(:) ! volumetric liquid water conetnt of snow at the start of the time step
+ real(dp), allocatable :: liqSoilInit(:) ! soil moisture at the start of the time step
+ ! sundials addition
+ logical(lgt) :: sundials=.false.
+ logical(lgt) :: tooMuchSublim ! flag to denote that there was too much sublimation in a given time step
+
+
+ ! ----------------------------------------------------------------------------------------------------------------------------------------------
+ ! initialize error control
+ err=0; message="coupled_em/"
! This is the start of a data step for a local HRU
- ! check that the decision is supported
- if(model_decisions(iLookDECISIONS%groundwatr)%iDecision==bigBucket .and. &
- model_decisions(iLookDECISIONS%spatial_gw)%iDecision/=localColumn)then
- message=trim(message)//'expect "spatial_gw" decision to equal localColumn when "groundwatr" decision is bigBucket'
- err=20; return
- endif
-
- ! check if the aquifer is included
- includeAquifer = (model_decisions(iLookDECISIONS%groundwatr)%iDecision==bigBucket)
-
- ! initialize the numerix tracking variables
- indx_data%var(iLookINDEX%numberFluxCalc )%dat(1) = 0 ! number of flux calculations (-)
- indx_data%var(iLookINDEX%numberStateSplit )%dat(1) = 0 ! number of state splitting solutions (-)
- indx_data%var(iLookINDEX%numberDomainSplitNrg )%dat(1) = 0 ! number of domain splitting solutions for energy (-)
- indx_data%var(iLookINDEX%numberDomainSplitMass)%dat(1) = 0 ! number of domain splitting solutions for mass (-)
- indx_data%var(iLookINDEX%numberScalarSolutions)%dat(1) = 0 ! number of scalar solutions (-)
-
- ! link canopy depth to the information in the data structure
- canopy: associate(canopyDepth => diag_data%var(iLookDIAG%scalarCanopyDepth)%dat(1) ) ! intent(out): [dp] canopy depth (m)
+ ! check that the decision is supported
+ if(model_decisions(iLookDECISIONS%groundwatr)%iDecision==bigBucket .and. &
+ model_decisions(iLookDECISIONS%spatial_gw)%iDecision/=localColumn)then
+ message=trim(message)//'expect "spatial_gw" decision to equal localColumn when "groundwatr" decision is bigBucket'
+ err=20; return
+ endif
+ ! check if the aquifer is included
+ includeAquifer = (model_decisions(iLookDECISIONS%groundwatr)%iDecision==bigBucket)
+
+ ! initialize the numerix tracking variables
+ indx_data%var(iLookINDEX%numberFluxCalc )%dat(1) = 0 ! number of flux calculations (-)
+ indx_data%var(iLookINDEX%numberStateSplit )%dat(1) = 0 ! number of state splitting solutions (-)
+ indx_data%var(iLookINDEX%numberDomainSplitNrg )%dat(1) = 0 ! number of domain splitting solutions for energy (-)
+ indx_data%var(iLookINDEX%numberDomainSplitMass)%dat(1) = 0 ! number of domain splitting solutions for mass (-)
+ indx_data%var(iLookINDEX%numberScalarSolutions)%dat(1) = 0 ! number of scalar solutions (-)
+
+ ! link canopy depth to the information in the data structure
+ canopy: associate(canopyDepth => diag_data%var(iLookDIAG%scalarCanopyDepth)%dat(1) ) ! intent(out): [dp] canopy depth (m)
+
- ! start by NOT pausing
- pauseFlag=.false.
+ ! start by NOT pausing
+ pauseFlag=.false.
- ! start by assuming that the step is successful
- stepFailure = .false.
- doLayerMerge = .false.
+ ! start by assuming that the step is successful
+ stepFailure = .false.
+ doLayerMerge = .false.
- ! initialize flags to modify the veg layers or modify snow layers
- modifiedLayers = .false. ! flag to denote that snow layers were modified
- modifiedVegState = .false. ! flag to denote that vegetation states were modified
+ ! initialize flags to modify the veg layers or modify snow layers
+ modifiedLayers = .false. ! flag to denote that snow layers were modified
+ modifiedVegState = .false. ! flag to denote that vegetation states were modified
- ! define the first step
- firstSubStep = .true.
+ ! define the first step
+ firstSubStep = .true.
- ! count the number of snow and soil layers
- ! NOTE: need to re-compute the number of snow and soil layers at the start of each sub-step because the number of layers may change
- ! (nSnow and nSoil are shared in the data structure)
- nSnow = count(indx_data%var(iLookINDEX%layerType)%dat==iname_snow)
- nSoil = count(indx_data%var(iLookINDEX%layerType)%dat==iname_soil)
+ ! count the number of snow and soil layers
+ ! NOTE: need to re-compute the number of snow and soil layers at the start of each sub-step because the number of layers may change
+ ! (nSnow and nSoil are shared in the data structure)
+ nSnow = count(indx_data%var(iLookINDEX%layerType)%dat==iname_snow)
+ nSoil = count(indx_data%var(iLookINDEX%layerType)%dat==iname_soil)
- ! compute the total number of snow and soil layers
- nLayers = nSnow + nSoil
-
+ ! compute the total number of snow and soil layers
+ nLayers = nSnow + nSoil
+
! create temporary data structures for prognostic variables
- call resizeData(prog_meta(:),prog_data,prog_temp,err=err,message=cmessage)
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
+ call resizeData(prog_meta(:),prog_data,prog_temp,err=err,message=cmessage)
+ if(err/=0)then
+ err=20
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ endif
- ! create temporary data structures for diagnostic variables
- call resizeData(diag_meta(:),diag_data,diag_temp,err=err,message=cmessage)
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
+ ! create temporary data structures for diagnostic variables
+ call resizeData(diag_meta(:),diag_data,diag_temp,err=err,message=cmessage)
+ if(err/=0)then
+ err=20
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ endif
- ! create temporary data structures for index variables
- call resizeData(indx_meta(:),indx_data,indx_temp,err=err,message=cmessage)
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
+ ! create temporary data structures for index variables
+ call resizeData(indx_meta(:),indx_data,indx_temp,err=err,message=cmessage)
+ if(err/=0)then
+ err=20
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ endif
- ! allocate space for the local fluxes
- call allocLocal(averageFlux_meta(:)%var_info,flux_mean,nSnow,nSoil,err,cmessage)
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if
+ ! allocate space for the local fluxes
+ call allocLocal(averageFlux_meta(:)%var_info,flux_mean,nSnow,nSoil,err,cmessage)
+ if(err/=0)then
+ err=20
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ end if
- ! initialize compression and surface melt pond
- sfcMeltPond = 0._dp ! change in storage associated with the surface melt pond (kg m-2)
- totalSoilCompress = 0._dp ! change in soil storage associated with compression of the matrix (kg m-2)
+ ! initialize compression and surface melt pond
+ sfcMeltPond = 0._dp ! change in storage associated with the surface melt pond (kg m-2)
+ totalSoilCompress = 0._dp ! change in soil storage associated with compression of the matrix (kg m-2)
- ! initialize mean fluxes
- do iVar=1,size(averageFlux_meta)
- flux_mean%var(iVar)%dat(:) = 0._dp
- end do
+ ! initialize mean fluxes
+ do iVar=1,size(averageFlux_meta)
+ flux_mean%var(iVar)%dat(:) = 0._dp
+ end do
- ! associate local variables with information in the data structures
- associate(&
- ! state variables in the vegetation canopy
- scalarCanopyLiq => prog_data%var(iLookPROG%scalarCanopyLiq)%dat(1) ,& ! canopy liquid water (kg m-2)
- scalarCanopyIce => prog_data%var(iLookPROG%scalarCanopyIce)%dat(1) ,& ! canopy ice content (kg m-2)
- ! state variables in the soil domain
- mLayerDepth => prog_data%var(iLookPROG%mLayerDepth)%dat(nSnow+1:nLayers) ,& ! depth of each soil layer (m)
- mLayerVolFracIce => prog_data%var(iLookPROG%mLayerVolFracIce)%dat(nSnow+1:nLayers) ,& ! volumetric ice content in each soil layer (-)
- mLayerVolFracLiq => prog_data%var(iLookPROG%mLayerVolFracLiq)%dat(nSnow+1:nLayers) ,& ! volumetric liquid water content in each soil layer (-)
- scalarAquiferStorage => prog_data%var(iLookPROG%scalarAquiferStorage)%dat(1) ,& ! aquifer storage (m)
- scalarTotalSoilIce => diag_data%var(iLookDIAG%scalarTotalSoilIce)%dat(1) ,& ! total ice in the soil column (kg m-2)
- scalarTotalSoilLiq => diag_data%var(iLookDIAG%scalarTotalSoilLiq)%dat(1) ,& ! total liquid water in the soil column (kg m-2)
- scalarTotalSoilWat => diag_data%var(iLookDIAG%scalarTotalSoilWat)%dat(1) & ! total water in the soil column (kg m-2)
- ) ! (association of local variables with information in the data structures
-
- ! save the liquid water and ice on the vegetation canopy
- scalarInitCanopyLiq = scalarCanopyLiq ! initial liquid water on the vegetation canopy (kg m-2)
- scalarInitCanopyIce = scalarCanopyIce ! initial ice on the vegetation canopy (kg m-2)
-
- ! compute total soil moisture and ice at the *START* of the step (kg m-2)
- scalarTotalSoilLiq = sum(iden_water*mLayerVolFracLiq(1:nSoil)*mLayerDepth(1:nSoil))
- scalarTotalSoilIce = sum(iden_water*mLayerVolFracIce(1:nSoil)*mLayerDepth(1:nSoil)) ! NOTE: no expansion and hence use iden_water
-
- ! compute storage of water in the canopy and the soil
- balanceCanopyWater0 = scalarCanopyLiq + scalarCanopyIce
- balanceSoilWater0 = scalarTotalSoilLiq + scalarTotalSoilIce
-
- ! get the total aquifer storage at the start of the time step (kg m-2)
- balanceAquifer0 = scalarAquiferStorage*iden_water
-
- ! save liquid water content
- if(printBalance)then
- allocate(liqSnowInit(nSnow), liqSoilInit(nSoil), stat=err)
- if(err/=0)then
- message=trim(message)//'unable to allocate space for the initial vectors'
- err=20; return
+ ! associate local variables with information in the data structures
+ associate(&
+ ! state variables in the vegetation canopy
+ scalarCanopyLiq => prog_data%var(iLookPROG%scalarCanopyLiq)%dat(1) ,& ! canopy liquid water (kg m-2)
+ scalarCanopyIce => prog_data%var(iLookPROG%scalarCanopyIce)%dat(1) ,& ! canopy ice content (kg m-2)
+ ! state variables in the soil domain
+ mLayerDepth => prog_data%var(iLookPROG%mLayerDepth)%dat(nSnow+1:nLayers) ,& ! depth of each soil layer (m)
+ mLayerVolFracIce => prog_data%var(iLookPROG%mLayerVolFracIce)%dat(nSnow+1:nLayers) ,& ! volumetric ice content in each soil layer (-)
+ mLayerVolFracLiq => prog_data%var(iLookPROG%mLayerVolFracLiq)%dat(nSnow+1:nLayers) ,& ! volumetric liquid water content in each soil layer (-)
+ scalarAquiferStorage => prog_data%var(iLookPROG%scalarAquiferStorage)%dat(1) ,& ! aquifer storage (m)
+ scalarTotalSoilIce => diag_data%var(iLookDIAG%scalarTotalSoilIce)%dat(1) ,& ! total ice in the soil column (kg m-2)
+ scalarTotalSoilLiq => diag_data%var(iLookDIAG%scalarTotalSoilLiq)%dat(1) ,& ! total liquid water in the soil column (kg m-2)
+ scalarTotalSoilWat => diag_data%var(iLookDIAG%scalarTotalSoilWat)%dat(1) & ! total water in the soil column (kg m-2)
+ ) ! (association of local variables with information in the data structures
+
+ ! save the liquid water and ice on the vegetation canopy
+ scalarInitCanopyLiq = scalarCanopyLiq ! initial liquid water on the vegetation canopy (kg m-2)
+ scalarInitCanopyIce = scalarCanopyIce ! initial ice on the vegetation canopy (kg m-2)
+
+ ! compute total soil moisture and ice at the *START* of the step (kg m-2)
+ scalarTotalSoilLiq = sum(iden_water*mLayerVolFracLiq(1:nSoil)*mLayerDepth(1:nSoil))
+ scalarTotalSoilIce = sum(iden_water*mLayerVolFracIce(1:nSoil)*mLayerDepth(1:nSoil)) ! NOTE: no expansion and hence use iden_water
+
+ ! compute storage of water in the canopy and the soil
+ balanceCanopyWater0 = scalarCanopyLiq + scalarCanopyIce
+ balanceSoilWater0 = scalarTotalSoilLiq + scalarTotalSoilIce
+
+ ! get the total aquifer storage at the start of the time step (kg m-2)
+ balanceAquifer0 = scalarAquiferStorage*iden_water
+
+ ! save liquid water content
+ if(printBalance)then
+ allocate(liqSnowInit(nSnow), liqSoilInit(nSoil), stat=err)
+ if(err/=0)then
+ message=trim(message)//'unable to allocate space for the initial vectors'
+ print*,message
+ err=20;
+ return
+ endif
+ if(nSnow>0) liqSnowInit = prog_data%var(iLookPROG%mLayerVolFracLiq)%dat(1:nSnow)
+ liqSoilInit = mLayerVolFracLiq
endif
- if(nSnow>0) liqSnowInit = prog_data%var(iLookPROG%mLayerVolFracLiq)%dat(1:nSnow)
- liqSoilInit = mLayerVolFracLiq
- endif
-
- ! end association of local variables with information in the data structures
- end associate
-
- ! short-cut to the algorithmic control parameters
- ! NOTE - temporary assignment of minstep to foce something reasonable
- minstep = 10._dp ! mpar_data%var(iLookPARAM%minstep)%dat(1) ! minimum time step (s)
- maxstep = mpar_data%var(iLookPARAM%maxstep)%dat(1) ! maximum time step (s)
- !print*, 'minstep, maxstep = ', minstep, maxstep
+ ! end association of local variables with information in the data structures
+ end associate
+
- ! compute the number of layers with roots
- nLayersRoots = count(prog_data%var(iLookPROG%iLayerHeight)%dat(nSnow:nLayers-1) < mpar_data%var(iLookPARAM%rootingDepth)%dat(1)-verySmall)
- if(nLayersRoots == 0)then
- message=trim(message)//'no roots within the soil profile'
- err=20; return
- end if
+ ! short-cut to the algorithmic control parameters
+ ! NOTE - temporary assignment of minstep to foce something reasonable
+ minstep = 10._dp ! mpar_data%var(iLookPARAM%minstep)%dat(1) ! minimum time step (s)
+ maxstep = mpar_data%var(iLookPARAM%maxstep)%dat(1) ! maximum time step (s)
+ !print*, 'minstep, maxstep = ', minstep, maxstep
+
+ ! compute the number of layers with roots
+ nLayersRoots = count(prog_data%var(iLookPROG%iLayerHeight)%dat(nSnow:nLayers-1) < mpar_data%var(iLookPARAM%rootingDepth)%dat(1)-verySmall)
+ if(nLayersRoots == 0)then
+ message=trim(message)//'no roots within the soil profile'
+ print*, message
+ err=20; return
+ end if
- ! define the foliage nitrogen factor
- diag_data%var(iLookDIAG%scalarFoliageNitrogenFactor)%dat(1) = 1._dp ! foliage nitrogen concentration (1.0 = saturated)
+ ! define the foliage nitrogen factor
+ diag_data%var(iLookDIAG%scalarFoliageNitrogenFactor)%dat(1) = 1._dp ! foliage nitrogen concentration (1.0 = saturated)
- ! save SWE
- oldSWE = prog_data%var(iLookPROG%scalarSWE)%dat(1)
- !print*, 'nSnow = ', nSnow
- !print*, 'oldSWE = ', oldSWE
+ ! save SWE
+ oldSWE = prog_data%var(iLookPROG%scalarSWE)%dat(1)
+ !print*, 'nSnow = ', nSnow
+ !print*, 'oldSWE = ', oldSWE
- ! *** compute phenology...
- ! ------------------------
-
+ ! *** compute phenology...
+ ! ------------------------
+
- ! compute the temperature of the root zone: used in vegetation phenology
- diag_data%var(iLookDIAG%scalarRootZoneTemp)%dat(1) = sum(prog_data%var(iLookPROG%mLayerTemp)%dat(nSnow+1:nSnow+nLayersRoots)) / real(nLayersRoots, kind(dp))
+ ! compute the temperature of the root zone: used in vegetation phenology
+ diag_data%var(iLookDIAG%scalarRootZoneTemp)%dat(1) = sum(prog_data%var(iLookPROG%mLayerTemp)%dat(nSnow+1:nSnow+nLayersRoots)) / real(nLayersRoots, kind(dp))
- ! remember if we compute the vegetation flux on the previous sub-step
- computeVegFluxOld = computeVegFlux
+ ! remember if we compute the vegetation flux on the previous sub-step
+ computeVegFluxOld = computeVegFlux
- ! compute the exposed LAI and SAI and whether veg is buried by snow
- call vegPhenlgy(&
- ! input/output: 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
- prog_data, & ! intent(in): model prognostic variables for a local HRU
- diag_data, & ! intent(inout): model diagnostic variables for a local HRU
- ! output
- computeVegFlux, & ! intent(out): flag to indicate if we are computing fluxes over vegetation (.false. means veg is buried with snow)
- canopyDepth, & ! intent(out): canopy depth (m)
- exposedVAI, & ! intent(out): exposed vegetation area index (m2 m-2)
- fracJulDay, & ! fractional julian days since the start of year
- yearLength, & ! number of days in the current year
- err,cmessage) ! intent(out): error control
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if
+ ! compute the exposed LAI and SAI and whether veg is buried by snow
+ call vegPhenlgy(&
+ ! input/output: 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
+ prog_data, & ! intent(in): model prognostic variables for a local HRU
+ diag_data, & ! intent(inout): model diagnostic variables for a local HRU
+ ! output
+ computeVegFlux, & ! intent(out): flag to indicate if we are computing fluxes over vegetation (.false. means veg is buried with snow)
+ canopyDepth, & ! intent(out): canopy depth (m)
+ exposedVAI, & ! intent(out): exposed vegetation area index (m2 m-2)
+ fracJulDay, & ! fractional julian days since the start of year
+ yearLength, & ! number of days in the current year
+ err,cmessage) ! intent(out): error control
+ if(err/=0)then
+ err=20
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ end if
- ! check
- if(computeVegFlux)then
- if(canopyDepth < epsilon(canopyDepth))then
- message=trim(message)//'canopy depth is zero when computeVegFlux flag is .true.'
- err=20; return
+ ! check
+ if(computeVegFlux)then
+ if(canopyDepth < epsilon(canopyDepth))then
+ message=trim(message)//'canopy depth is zero when computeVegFlux flag is .true.'
+ print*, message
+ err=20; return
+ endif
endif
- endif
- ! flag the case where number of vegetation states has changed
- modifiedVegState = (computeVegFlux.neqv.computeVegFluxOld)
-
- ! *** compute wetted canopy area...
- ! ---------------------------------
-
- ! compute maximum canopy liquid water (kg m-2)
- diag_data%var(iLookDIAG%scalarCanopyLiqMax)%dat(1) = mpar_data%var(iLookPARAM%refInterceptCapRain)%dat(1)*exposedVAI
-
- ! compute maximum canopy ice content (kg m-2)
- ! NOTE 1: this is used to compute the snow fraction on the canopy, as used in *BOTH* the radiation AND canopy sublimation routines
- ! NOTE 2: this is a different variable than the max ice used in the throughfall (snow interception) calculations
- ! NOTE 3: use maximum per unit leaf area storage capacity for snow (kg m-2)
- select case(model_decisions(iLookDECISIONS%snowIncept)%iDecision)
- case(lightSnow); diag_data%var(iLookDIAG%scalarCanopyIceMax)%dat(1) = exposedVAI*mpar_data%var(iLookPARAM%refInterceptCapSnow)%dat(1)
- case(stickySnow); diag_data%var(iLookDIAG%scalarCanopyIceMax)%dat(1) = exposedVAI*mpar_data%var(iLookPARAM%refInterceptCapSnow)%dat(1)*4._dp
- case default; message=trim(message)//'unable to identify option for maximum branch interception capacity'; err=20; return
- end select ! identifying option for maximum branch interception capacity
- !print*, 'diag_data%var(iLookDIAG%scalarCanopyLiqMax)%dat(1) = ', diag_data%var(iLookDIAG%scalarCanopyLiqMax)%dat(1)
- !print*, 'diag_data%var(iLookDIAG%scalarCanopyIceMax)%dat(1) = ', diag_data%var(iLookDIAG%scalarCanopyIceMax)%dat(1)
-
- ! compute wetted fraction of the canopy
- ! NOTE: assume that the wetted fraction is constant over the substep for the radiation calculations
- if(computeVegFlux)then
+ ! flag the case where number of vegetation states has changed
+ modifiedVegState = (computeVegFlux.neqv.computeVegFluxOld)
+
+ ! *** compute wetted canopy area...
+ ! ---------------------------------
+
+ ! compute maximum canopy liquid water (kg m-2)
+ diag_data%var(iLookDIAG%scalarCanopyLiqMax)%dat(1) = mpar_data%var(iLookPARAM%refInterceptCapRain)%dat(1)*exposedVAI
+
+ ! compute maximum canopy ice content (kg m-2)
+ ! NOTE 1: this is used to compute the snow fraction on the canopy, as used in *BOTH* the radiation AND canopy sublimation routines
+ ! NOTE 2: this is a different variable than the max ice used in the throughfall (snow interception) calculations
+ ! NOTE 3: use maximum per unit leaf area storage capacity for snow (kg m-2)
+ select case(model_decisions(iLookDECISIONS%snowIncept)%iDecision)
+ case(lightSnow); diag_data%var(iLookDIAG%scalarCanopyIceMax)%dat(1) = exposedVAI*mpar_data%var(iLookPARAM%refInterceptCapSnow)%dat(1)
+ case(stickySnow); diag_data%var(iLookDIAG%scalarCanopyIceMax)%dat(1) = exposedVAI*mpar_data%var(iLookPARAM%refInterceptCapSnow)%dat(1)*4._dp
+ case default
+ message=trim(message)//'unable to identify option for maximum branch interception capacity'
+ print*, message
+ err=20
+ return
+ end select ! identifying option for maximum branch interception capacity
+ !print*, 'diag_data%var(iLookDIAG%scalarCanopyLiqMax)%dat(1) = ', diag_data%var(iLookDIAG%scalarCanopyLiqMax)%dat(1)
+ !print*, 'diag_data%var(iLookDIAG%scalarCanopyIceMax)%dat(1) = ', diag_data%var(iLookDIAG%scalarCanopyIceMax)%dat(1)
! compute wetted fraction of the canopy
- call wettedFrac(&
+ ! NOTE: assume that the wetted fraction is constant over the substep for the radiation calculations
+ if(computeVegFlux)then
+
+ ! compute wetted fraction of the canopy
+ call wettedFrac(&
! input
.false., & ! flag to denote if derivatives are required
.false., & ! flag to denote if derivatives are calculated numerically
@@ -464,73 +502,92 @@ subroutine coupled_em(&
dCanopyWetFraction_dWat, & ! derivative in wetted fraction w.r.t. canopy liquid water content (kg-1 m2)
dCanopyWetFraction_dT, & ! derivative in wetted fraction w.r.t. canopy liquid water content (kg-1 m2)
err,cmessage)
- if(err/=0)then; message=trim(message)//trim(cmessage); return; end if
+ if(err/=0)then
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ end if
- ! vegetation is completely buried by snow (or no veg exists at all)
- else
- diag_data%var(iLookDIAG%scalarCanopyWetFraction)%dat(1) = 0._dp
- dCanopyWetFraction_dWat = 0._dp
- dCanopyWetFraction_dT = 0._dp
- end if
-
+ ! vegetation is completely buried by snow (or no veg exists at all)
+ else
+ diag_data%var(iLookDIAG%scalarCanopyWetFraction)%dat(1) = 0._dp
+ dCanopyWetFraction_dWat = 0._dp
+ dCanopyWetFraction_dT = 0._dp
+ end if
+
- ! *** compute snow albedo...
- ! --------------------------
- ! NOTE: this should be done before the radiation calculations
- ! NOTE: uses snowfall; should really use canopy throughfall + canopy unloading
- call snowAlbedo(&
- ! input: model control
- data_step, & ! intent(in): model time step (s)
- (nSnow > 0), & ! intent(in): logical flag to denote if snow is present
- ! input/output: data structures
- model_decisions, & ! intent(in): model decisions
- mpar_data, & ! intent(in): model parameters
- flux_data, & ! intent(in): model flux variables
- diag_data, & ! intent(inout): model diagnostic variables for a local HRU
- prog_data, & ! intent(inout): model prognostic variables for a local HRU
- ! output: error control
- err,cmessage) ! intent(out): error control
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if
+ ! *** compute snow albedo...
+ ! --------------------------
+ ! NOTE: this should be done before the radiation calculations
+ ! NOTE: uses snowfall; should really use canopy throughfall + canopy unloading
+ call snowAlbedo(&
+ ! input: model control
+ data_step, & ! intent(in): model time step (s)
+ (nSnow > 0), & ! intent(in): logical flag to denote if snow is present
+ ! input/output: data structures
+ model_decisions, & ! intent(in): model decisions
+ mpar_data, & ! intent(in): model parameters
+ flux_data, & ! intent(in): model flux variables
+ diag_data, & ! intent(inout): model diagnostic variables for a local HRU
+ prog_data, & ! intent(inout): model prognostic variables for a local HRU
+ ! output: error control
+ err,cmessage) ! intent(out): error control
+ if(err/=0)then
+ err=20
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ end if
- ! *** compute canopy sw radiation fluxes...
- ! -----------------------------------------
- call vegSWavRad(&
- data_step, & ! intent(in): time step (s) -- only used in Noah-MP radiation, to compute albedo
- nSnow, & ! intent(in): number of snow layers
- nSoil, & ! intent(in): number of soil layers
- nLayers, & ! intent(in): total number of layers
- computeVegFlux, & ! intent(in): logical flag to compute vegetation fluxes (.false. if veg buried by snow)
- type_data, & ! intent(in): type of vegetation and soil
- prog_data, & ! intent(inout): model prognostic variables for a local HRU
- diag_data, & ! intent(inout): model diagnostic variables for a local HRU
- flux_data, & ! intent(inout): model flux variables
- err,cmessage) ! intent(out): error control
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if
+ ! *** compute canopy sw radiation fluxes...
+ ! -----------------------------------------
+ call vegSWavRad(&
+ data_step, & ! intent(in): time step (s) -- only used in Noah-MP radiation, to compute albedo
+ nSnow, & ! intent(in): number of snow layers
+ nSoil, & ! intent(in): number of soil layers
+ nLayers, & ! intent(in): total number of layers
+ computeVegFlux, & ! intent(in): logical flag to compute vegetation fluxes (.false. if veg buried by snow)
+ type_data, & ! intent(in): type of vegetation and soil
+ prog_data, & ! intent(inout): model prognostic variables for a local HRU
+ diag_data, & ! intent(inout): model diagnostic variables for a local HRU
+ flux_data, & ! intent(inout): model flux variables
+ err,cmessage) ! intent(out): error control
+ if(err/=0)then
+ err=20
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ end if
- ! *** compute canopy throughfall and unloading...
- ! -----------------------------------------------
- ! NOTE 1: this needs to be done before solving the energy and liquid water equations, to account for the heat advected with precipitation (and throughfall/unloading)
- ! NOTE 2: the unloading flux is computed using canopy drip (scalarCanopyLiqDrainage) from the previous time step
- call canopySnow(&
- ! input: model control
- data_step, & ! intent(in): time step (seconds)
- exposedVAI, & ! intent(in): exposed vegetation area index (m2 m-2)
- computeVegFlux, & ! intent(in): flag to denote if computing energy flux over vegetation
- ! input/output: data structures
- model_decisions, & ! intent(in): model decisions
- forc_data, & ! intent(in): model forcing data
- mpar_data, & ! intent(in): model parameters
- diag_data, & ! intent(in): model diagnostic variables for a local HRU
- prog_data, & ! intent(inout): model prognostic variables for a local HRU
- flux_data, & ! intent(inout): model flux variables
- ! output: error control
- err,cmessage) ! intent(out): error control
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if
+ ! *** compute canopy throughfall and unloading...
+ ! -----------------------------------------------
+ ! NOTE 1: this needs to be done before solving the energy and liquid water equations, to account for the heat advected with precipitation (and throughfall/unloading)
+ ! NOTE 2: the unloading flux is computed using canopy drip (scalarCanopyLiqDrainage) from the previous time step
+ call canopySnow(&
+ ! input: model control
+ data_step, & ! intent(in): time step (seconds)
+ exposedVAI, & ! intent(in): exposed vegetation area index (m2 m-2)
+ computeVegFlux, & ! intent(in): flag to denote if computing energy flux over vegetation
+ ! input/output: data structures
+ model_decisions, & ! intent(in): model decisions
+ forc_data, & ! intent(in): model forcing data
+ mpar_data, & ! intent(in): model parameters
+ diag_data, & ! intent(in): model diagnostic variables for a local HRU
+ prog_data, & ! intent(inout): model prognostic variables for a local HRU
+ flux_data, & ! intent(inout): model flux variables
+ ! output: error control
+ err,cmessage) ! intent(out): error control
+ if(err/=0)then
+ err=20
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ end if
- ! adjust canopy temperature to account for new snow
- if(computeVegFlux)then ! logical flag to compute vegetation fluxes (.false. if veg buried by snow)
- call tempAdjust(&
+ ! adjust canopy temperature to account for new snow
+ if(computeVegFlux)then ! logical flag to compute vegetation fluxes (.false. if veg buried by snow)
+ call tempAdjust(&
! input: derived parameters
canopyDepth, & ! intent(in): canopy depth (m)
! input/output: data structures
@@ -539,161 +596,199 @@ subroutine coupled_em(&
diag_data, & ! intent(out): model diagnostic variables for a local HRU
! output: error control
err,cmessage) ! intent(out): error control
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if
+ if(err/=0)then
+ err=20
+ message=trim(message)//trim(cmessage)
+ return
+ end if
endif ! if computing fluxes over vegetation
- ! initialize drainage and throughfall
- ! NOTE 1: this needs to be done before solving the energy and liquid water equations, to account for the heat advected with precipitation
- ! NOTE 2: this initialization needs to be done AFTER the call to canopySnow, since canopySnow uses canopy drip drom the previous time step
- if(.not.computeVegFlux)then
- flux_data%var(iLookFLUX%scalarThroughfallRain)%dat(1) = flux_data%var(iLookFLUX%scalarRainfall)%dat(1)
- flux_data%var(iLookFLUX%scalarCanopyLiqDrainage)%dat(1) = 0._dp
- else
- flux_data%var(iLookFLUX%scalarThroughfallRain)%dat(1) = 0._dp
- flux_data%var(iLookFLUX%scalarCanopyLiqDrainage)%dat(1) = 0._dp
- end if
-
- ! ****************************************************************************************************
- ! *** MAIN SOLVER ************************************************************************************
- ! ****************************************************************************************************
-
- ! initialize the length of the sub-step
- dt_solv = 0._dp ! length of time step that has been completed (s)
- dt_init = min(data_step,maxstep) / dt_init_factor ! initial substep length (s)
- dt_sub = dt_init ! length of substep
- dtSave = dt_init ! length of substep
-
- ! initialize the number of sub-steps
- nsub=0
-
- ! loop through sub-steps
- substeps: do ! continuous do statement with exit clause (alternative to "while")
-
- ! print progress
- !print*, '*** new substep'
- !write(*,'(a,3(f11.4,1x))') 'dt_sub, dt_init = ', dt_sub, dt_init
-
- ! print progress
- if(globalPrintFlag)then
- write(*,'(a,1x,4(f13.5,1x))') ' start of step: dt_init, dt_sub, dt_solv, data_step: ', dt_init, dt_sub, dt_solv, data_step
- print*, 'stepFailure = ', stepFailure
- print*, 'before resizeData: nSnow, nSoil = ', nSnow, nSoil
- endif
-
- ! increment the number of sub-steps
- nsub = nsub+1
-
- ! resize the "indx_data" structure
- ! NOTE: this is necessary because the length of index variables depends on a given split
- ! --> the resize here is overwritten later (in indexSplit)
- ! --> admittedly ugly, and retained for now
- if(stepFailure)then
- call resizeData(indx_meta(:),indx_temp,indx_data,err=err,message=cmessage)
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
+ ! initialize drainage and throughfall
+ ! NOTE 1: this needs to be done before solving the energy and liquid water equations, to account for the heat advected with precipitation
+ ! NOTE 2: this initialization needs to be done AFTER the call to canopySnow, since canopySnow uses canopy drip drom the previous time step
+ if(.not.computeVegFlux)then
+ flux_data%var(iLookFLUX%scalarThroughfallRain)%dat(1) = flux_data%var(iLookFLUX%scalarRainfall)%dat(1)
+ flux_data%var(iLookFLUX%scalarCanopyLiqDrainage)%dat(1) = 0._dp
else
- call resizeData(indx_meta(:),indx_data,indx_temp,err=err,message=cmessage)
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
- endif
-
- ! save/recover copies of index variables
- do iVar=1,size(indx_data%var)
- !print*, 'indx_meta(iVar)%varname = ', trim(indx_meta(iVar)%varname)
- select case(stepFailure)
- case(.false.); indx_temp%var(iVar)%dat(:) = indx_data%var(iVar)%dat(:)
- case(.true.); indx_data%var(iVar)%dat(:) = indx_temp%var(iVar)%dat(:)
- end select
- end do ! looping through variables
-
- ! save/recover copies of prognostic variables
- do iVar=1,size(prog_data%var)
- !print*, 'prog_meta(iVar)%varname = ', trim(prog_meta(iVar)%varname)
- select case(stepFailure)
- case(.false.); prog_temp%var(iVar)%dat(:) = prog_data%var(iVar)%dat(:)
- case(.true.); prog_data%var(iVar)%dat(:) = prog_temp%var(iVar)%dat(:)
- end select
- end do ! looping through variables
-
- ! save/recover copies of diagnostic variables
- do iVar=1,size(diag_data%var)
- select case(stepFailure)
- case(.false.); diag_temp%var(iVar)%dat(:) = diag_data%var(iVar)%dat(:)
- case(.true.); diag_data%var(iVar)%dat(:) = diag_temp%var(iVar)%dat(:)
- end select
- end do ! looping through variables
-
- ! re-assign dimension lengths
- nSnow = count(indx_data%var(iLookINDEX%layerType)%dat==iname_snow)
- nSoil = count(indx_data%var(iLookINDEX%layerType)%dat==iname_soil)
- nLayers = nSnow+nSoil
-
-
- ! *** merge/sub-divide snow layers...
- ! -----------------------------------
- call volicePack(&
- ! input/output: model data structures
- doLayerMerge, & ! intent(in): flag to force merge of snow layers
- model_decisions, & ! intent(in): model decisions
- mpar_data, & ! intent(in): model parameters
- indx_data, & ! intent(inout): type of each layer
- prog_data, & ! intent(inout): model prognostic variables for a local HRU
- diag_data, & ! intent(inout): model diagnostic variables for a local HRU
- flux_data, & ! intent(inout): model fluxes for a local HRU
- ! output
- modifiedLayers, & ! intent(out): flag to denote that layers were modified
- err,cmessage) ! intent(out): error control
- if(err/=0)then; err=55; message=trim(message)//trim(cmessage); return; end if
-
- ! save the number of snow and soil layers
- nSnow = indx_data%var(iLookINDEX%nSnow)%dat(1)
- nSoil = indx_data%var(iLookINDEX%nSoil)%dat(1)
- nLayers = indx_data%var(iLookINDEX%nLayers)%dat(1)
-
-
- ! compute the indices for the model state variables
- if(firstSubStep .or. modifiedVegState .or. modifiedLayers)then
- call indexState(computeVegFlux, & ! intent(in): flag to denote if computing the vegetation flux
- includeAquifer, & ! intent(in): flag to denote if included the aquifer
- nSnow,nSoil,nLayers, & ! intent(in): number of snow and soil layers, and total number of layers
- indx_data, & ! intent(inout): indices defining model states and layers
- err,cmessage) ! intent(out): error control
- if(err/=0)then; message=trim(message)//trim(cmessage); return; end if
+ flux_data%var(iLookFLUX%scalarThroughfallRain)%dat(1) = 0._dp
+ flux_data%var(iLookFLUX%scalarCanopyLiqDrainage)%dat(1) = 0._dp
end if
- ! recreate the temporary data structures
- ! NOTE: resizeData(meta, old, new, ..)
- if(modifiedVegState .or. modifiedLayers)then
-
- ! create temporary data structures for prognostic variables
- call resizeData(prog_meta(:),prog_data,prog_temp,copy=.true.,err=err,message=cmessage)
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
-
- ! create temporary data structures for diagnostic variables
- call resizeData(diag_meta(:),diag_data,diag_temp,copy=.true.,err=err,message=cmessage)
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
-
- ! create temporary data structures for index variables
- call resizeData(indx_meta(:),indx_data,indx_temp,copy=.true.,err=err,message=cmessage)
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
-
- do iVar=1,size(indx_data%var)
- !print*, 'indx_meta(iVar)%varname = ', trim(indx_meta(iVar)%varname)
- select case(stepFailure)
- case(.false.); indx_temp%var(iVar)%dat(:) = indx_data%var(iVar)%dat(:)
- case(.true.); indx_data%var(iVar)%dat(:) = indx_temp%var(iVar)%dat(:)
- end select
- end do ! looping through variables
+ ! ****************************************************************************************************
+ ! *** MAIN SOLVER ************************************************************************************
+ ! ****************************************************************************************************
- endif ! if modified the states
+ ! initialize the length of the sub-step
+ dt_solv = 0._dp ! length of time step that has been completed (s)
+ dt_init = min(data_step,maxstep) / dt_init_factor ! initial substep length (s)
+ dt_sub = dt_init ! length of substep
+ dtSave = dt_init ! length of substep
+
+ ! initialize the number of sub-steps
+ nsub=0
+
+ ! loop through sub-steps
+ substeps: do ! continuous do statement with exit clause (alternative to "while")
+
+ ! print progress
+ !print*, '*** new substep'
+ !write(*,'(a,3(f11.4,1x))') 'dt_sub, dt_init = ', dt_sub, dt_init
+
+ ! print progress
+ if(globalPrintFlag)then
+ write(*,'(a,1x,4(f13.5,1x))') ' start of step: dt_init, dt_sub, dt_solv, data_step: ', dt_init, dt_sub, dt_solv, data_step
+ print*, 'stepFailure = ', stepFailure
+ print*, 'before resizeData: nSnow, nSoil = ', nSnow, nSoil
+ endif
+
+ ! increment the number of sub-steps
+ nsub = nsub+1
+
+ ! resize the "indx_data" structure
+ ! NOTE: this is necessary because the length of index variables depends on a given split
+ ! --> the resize here is overwritten later (in indexSplit)
+ ! --> admittedly ugly, and retained for now
+ if(stepFailure)then
+ call resizeData(indx_meta(:),indx_temp,indx_data,err=err,message=cmessage)
+ if(err/=0)then
+ err=20
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ endif
+ else
+ call resizeData(indx_meta(:),indx_data,indx_temp,err=err,message=cmessage)
+ if(err/=0)then
+ err=20
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ endif
+ endif
+
+ ! save/recover copies of index variables
+ do iVar=1,size(indx_data%var)
+ !print*, 'indx_meta(iVar)%varname = ', trim(indx_meta(iVar)%varname)
+ select case(stepFailure)
+ case(.false.); indx_temp%var(iVar)%dat(:) = indx_data%var(iVar)%dat(:)
+ case(.true.); indx_data%var(iVar)%dat(:) = indx_temp%var(iVar)%dat(:)
+ end select
+ end do ! looping through variables
+
+ ! save/recover copies of prognostic variables
+ do iVar=1,size(prog_data%var)
+ !print*, 'prog_meta(iVar)%varname = ', trim(prog_meta(iVar)%varname)
+ select case(stepFailure)
+ case(.false.); prog_temp%var(iVar)%dat(:) = prog_data%var(iVar)%dat(:)
+ case(.true.); prog_data%var(iVar)%dat(:) = prog_temp%var(iVar)%dat(:)
+ end select
+ end do ! looping through variables
+
+ ! save/recover copies of diagnostic variables
+ do iVar=1,size(diag_data%var)
+ select case(stepFailure)
+ case(.false.); diag_temp%var(iVar)%dat(:) = diag_data%var(iVar)%dat(:)
+ case(.true.); diag_data%var(iVar)%dat(:) = diag_temp%var(iVar)%dat(:)
+ end select
+ end do ! looping through variables
+
+ ! re-assign dimension lengths
+ nSnow = count(indx_data%var(iLookINDEX%layerType)%dat==iname_snow)
+ nSoil = count(indx_data%var(iLookINDEX%layerType)%dat==iname_soil)
+ nLayers = nSnow+nSoil
+
- ! define the number of state variables
- nState = indx_data%var(iLookINDEX%nState)%dat(1)
+ ! *** merge/sub-divide snow layers...
+ ! -----------------------------------
+ call volicePack(&
+ ! input/output: model data structures
+ doLayerMerge, & ! intent(in): flag to force merge of snow layers
+ model_decisions, & ! intent(in): model decisions
+ mpar_data, & ! intent(in): model parameters
+ indx_data, & ! intent(inout): type of each layer
+ prog_data, & ! intent(inout): model prognostic variables for a local HRU
+ diag_data, & ! intent(inout): model diagnostic variables for a local HRU
+ flux_data, & ! intent(inout): model fluxes for a local HRU
+ ! output
+ modifiedLayers, & ! intent(out): flag to denote that layers were modified
+ err,cmessage) ! intent(out): error control
+ if(err/=0)then
+ err=55
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ end if
+
+ ! save the number of snow and soil layers
+ nSnow = indx_data%var(iLookINDEX%nSnow)%dat(1)
+ nSoil = indx_data%var(iLookINDEX%nSoil)%dat(1)
+ nLayers = indx_data%var(iLookINDEX%nLayers)%dat(1)
+
+
+ ! compute the indices for the model state variables
+ if(firstSubStep .or. modifiedVegState .or. modifiedLayers)then
+ call indexState(computeVegFlux, & ! intent(in): flag to denote if computing the vegetation flux
+ includeAquifer, & ! intent(in): flag to denote if included the aquifer
+ nSnow,nSoil,nLayers, & ! intent(in): number of snow and soil layers, and total number of layers
+ indx_data, & ! intent(inout): indices defining model states and layers
+ err,cmessage) ! intent(out): error control
+ if(err/=0)then
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ end if
+ end if
+
+ ! recreate the temporary data structures
+ ! NOTE: resizeData(meta, old, new, ..)
+ if(modifiedVegState .or. modifiedLayers)then
+
+ ! create temporary data structures for prognostic variables
+ call resizeData(prog_meta(:),prog_data,prog_temp,copy=.true.,err=err,message=cmessage)
+ if(err/=0)then;
+ err=20
+ message=trim(message)//trim(cmessage);
+ print*, message
+ return
+ endif
+
+ ! create temporary data structures for diagnostic variables
+ call resizeData(diag_meta(:),diag_data,diag_temp,copy=.true.,err=err,message=cmessage)
+ if(err/=0)then
+ err=20
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ endif
+
+ ! create temporary data structures for index variables
+ call resizeData(indx_meta(:),indx_data,indx_temp,copy=.true.,err=err,message=cmessage)
+ if(err/=0)then
+ err=20
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ endif
+
+ do iVar=1,size(indx_data%var)
+ !print*, 'indx_meta(iVar)%varname = ', trim(indx_meta(iVar)%varname)
+ select case(stepFailure)
+ case(.false.); indx_temp%var(iVar)%dat(:) = indx_data%var(iVar)%dat(:)
+ case(.true.); indx_data%var(iVar)%dat(:) = indx_temp%var(iVar)%dat(:)
+ end select
+ end do ! looping through variables
+
+ endif ! if modified the states
+
+ ! define the number of state variables
+ nState = indx_data%var(iLookINDEX%nState)%dat(1)
- ! *** compute diagnostic variables for each layer...
- ! --------------------------------------------------
- ! NOTE: this needs to be done AFTER volicePack, since layers may have been sub-divided and/or merged
- call diagn_evar(&
+ ! *** compute diagnostic variables for each layer...
+ ! --------------------------------------------------
+ ! NOTE: this needs to be done AFTER volicePack, since layers may have been sub-divided and/or merged
+ call diagn_evar(&
! input: control variables
computeVegFlux, & ! intent(in): flag to denote if computing the vegetation flux
canopyDepth, & ! intent(in): canopy depth (m)
@@ -704,15 +799,19 @@ subroutine coupled_em(&
diag_data, & ! intent(inout): model diagnostic variables for a local HRU
! output: error control
err,cmessage) ! intent(out): error control
- if(err/=0)then; err=55; message=trim(message)//trim(cmessage); return; end if
-
-
- ! *** compute melt of the "snow without a layer"...
- ! -------------------------------------------------
- ! NOTE: forms a surface melt pond, which drains into the upper-most soil layer through the time step
- ! (check for the special case of "snow without a layer")
- if(nSnow==0)then
- call implctMelt(&
+ if(err/=0)then
+ err=55;
+ message=trim(message)//trim(cmessage)
+ return
+ end if
+
+
+ ! *** compute melt of the "snow without a layer"...
+ ! -------------------------------------------------
+ ! NOTE: forms a surface melt pond, which drains into the upper-most soil layer through the time step
+ ! (check for the special case of "snow without a layer")
+ if(nSnow==0)then
+ call implctMelt(&
! input/output: integrated snowpack properties
prog_data%var(iLookPROG%scalarSWE)%dat(1), & ! intent(inout): snow water equivalent (kg m-2)
prog_data%var(iLookPROG%scalarSnowDepth)%dat(1), & ! intent(inout): snow depth (m)
@@ -723,8 +822,13 @@ subroutine coupled_em(&
diag_data%var(iLookDIAG%mLayerVolHtCapBulk)%dat(nSnow+1),& ! intent(inout): surface layer volumetric heat capacity (J m-3 K-1)
! output: error control
err,cmessage ) ! intent(out): error control
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if
- end if
+ if(err/=0)then
+ err=20
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ end if
+ end if ! nsnow == 0
! *** solve model equations...
! ----------------------------
@@ -762,34 +866,40 @@ subroutine coupled_em(&
err,cmessage) ! intent(out): error code and error message
! check for all errors (error recovery within opSplittin)
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if
+ if(err/=0)then
+ err=20
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ end if
! process the flag for too much melt
if(tooMuchMelt)then
- stepFailure = .true.
- doLayerMerge = .true.
+ stepFailure = .true.
+ doLayerMerge = .true.
else
- doLayerMerge = .false.
+ doLayerMerge = .false.
endif
! handle special case of the step failure
! NOTE: need to revert back to the previous state vector that we were happy with and reduce the time step
+ ! TODO: ask isn't this what the actors program does without the code block below
if(stepFailure)then
- ! halve step
- dt_sub = dtSave/2._dp
+ ! halve step
+ dt_sub = dtSave/2._dp
- ! check that the step is not tiny
- if(dt_sub < minstep)then
- print*,ixSolution
- print*, 'dtSave, dt_sub', dtSave, dt_sub
- message=trim(message)//'length of the coupled step is below the minimum step length'
- err=20; return
- endif
+ ! check that the step is not tiny
+ if(dt_sub < minstep)then
+ print*,ixSolution
+ print*, 'dtSave, dt_sub', dtSave, dt_sub
+ message=trim(message)//'length of the coupled step is below the minimum step length'
+ err=20; return
+ endif
- ! try again
- cycle substeps
+ ! try again
+ cycle substeps
endif
@@ -816,124 +926,224 @@ subroutine coupled_em(&
! ------------------------------------------------------------
if(computeVegFlux)then
- ! remove mass of ice on the canopy
- scalarCanopyIce = scalarCanopyIce + scalarCanopySublimation*dt_sub
-
- ! if removed all ice, take the remaining sublimation from water
- if(scalarCanopyIce < 0._dp)then
- scalarCanopyLiq = scalarCanopyLiq + scalarCanopyIce
- scalarCanopyIce = 0._dp
- endif
-
- ! modify fluxes if there is insufficient canopy water to support the converged sublimation rate over the time step dt_sub
- if(scalarCanopyLiq < 0._dp)then
- ! --> superfluous sublimation flux
- superflousSub = -scalarCanopyLiq/dt_sub ! kg m-2 s-1
- superflousNrg = superflousSub*LH_sub ! W m-2 (J m-2 s-1)
- ! --> update fluxes and states
- scalarCanopySublimation = scalarCanopySublimation + superflousSub
- scalarLatHeatCanopyEvap = scalarLatHeatCanopyEvap + superflousNrg
- scalarSenHeatCanopy = scalarSenHeatCanopy - superflousNrg
- scalarCanopyLiq = 0._dp
- endif
+ ! remove mass of ice on the canopy
+ scalarCanopyIce = scalarCanopyIce + scalarCanopySublimation*dt_sub
+
+ ! if removed all ice, take the remaining sublimation from water
+ if(scalarCanopyIce < 0._dp)then
+ scalarCanopyLiq = scalarCanopyLiq + scalarCanopyIce
+ scalarCanopyIce = 0._dp
+ endif
+
+ ! modify fluxes if there is insufficient canopy water to support the converged sublimation rate over the time step dt_sub
+ if(scalarCanopyLiq < 0._dp)then
+ ! --> superfluous sublimation flux
+ superflousSub = -scalarCanopyLiq/dt_sub ! kg m-2 s-1
+ superflousNrg = superflousSub*LH_sub ! W m-2 (J m-2 s-1)
+ ! --> update fluxes and states
+ scalarCanopySublimation = scalarCanopySublimation + superflousSub
+ scalarLatHeatCanopyEvap = scalarLatHeatCanopyEvap + superflousNrg
+ scalarSenHeatCanopy = scalarSenHeatCanopy - superflousNrg
+ scalarCanopyLiq = 0._dp
+ endif
end if ! (if computing the vegetation flux)
- ! * compute change in ice content of the top snow layer due to sublimation...
- ! ---------------------------------------------------------------------------
- ! NOTE: this is done BEFORE densification
- if(nSnow > 0)then ! snow layers exist
-
- ! try to remove ice from the top layer
- iSnow=1
-
- ! save the mass of liquid water (kg m-2)
- massLiquid = mLayerDepth(iSnow)*mLayerVolFracLiq(iSnow)*iden_water
-
- ! add/remove the depth of snow gained/lost by frost/sublimation (m)
- ! NOTE: assume constant density
- mLayerDepth(iSnow) = mLayerDepth(iSnow) + dt_sub*scalarSnowSublimation/(mLayerVolFracIce(iSnow)*iden_ice)
-
- ! check that we did not remove the entire layer
- if(mLayerDepth(iSnow) < verySmall)then
- stepFailure = .true.
- doLayerMerge = .true.
- dt_sub = max(dtSave/2._dp, minstep)
- cycle substeps
- else
- stepFailure = .false.
- doLayerMerge = .false.
- endif
-
- ! update the volumetric fraction of liquid water
- mLayerVolFracLiq(iSnow) = massLiquid / (mLayerDepth(iSnow)*iden_water)
-
- ! no snow
- else
-
- ! no snow: check that sublimation is zero
- if(abs(scalarSnowSublimation) > verySmall)then
- message=trim(message)//'sublimation of snow has been computed when no snow exists'
- err=20; return
- end if
-
- end if ! (if snow layers exist)
-
- end associate sublime
-
- ! *** account for compaction and cavitation in the snowpack...
- ! ------------------------------------------------------------
- if(nSnow>0)then
- call snwDensify(&
- ! intent(in): variables
- dt_sub, & ! intent(in): time step (s)
- indx_data%var(iLookINDEX%nSnow)%dat(1), & ! intent(in): number of snow layers
- prog_data%var(iLookPROG%mLayerTemp)%dat(1:nSnow), & ! intent(in): temperature of each layer (K)
- diag_data%var(iLookDIAG%mLayerMeltFreeze)%dat(1:nSnow), & ! intent(in): volumetric melt in each layer (kg m-3)
- ! intent(in): parameters
- mpar_data%var(iLookPARAM%densScalGrowth)%dat(1), & ! intent(in): density scaling factor for grain growth (kg-1 m3)
- mpar_data%var(iLookPARAM%tempScalGrowth)%dat(1), & ! intent(in): temperature scaling factor for grain growth (K-1)
- mpar_data%var(iLookPARAM%grainGrowthRate)%dat(1), & ! intent(in): rate of grain growth (s-1)
- mpar_data%var(iLookPARAM%densScalOvrbdn)%dat(1), & ! intent(in): density scaling factor for overburden pressure (kg-1 m3)
- mpar_data%var(iLookPARAM%tempScalOvrbdn)%dat(1), & ! intent(in): temperature scaling factor for overburden pressure (K-1)
- mpar_data%var(iLookPARAM%baseViscosity)%dat(1), & ! intent(in): viscosity coefficient at T=T_frz and snow density=0 (kg m-2 s)
- ! intent(inout): state variables
- prog_data%var(iLookPROG%mLayerDepth)%dat(1:nSnow), & ! intent(inout): depth of each layer (m)
- prog_data%var(iLookPROG%mLayerVolFracLiq)%dat(1:nSnow), & ! intent(inout): volumetric fraction of liquid water after itertations (-)
- prog_data%var(iLookPROG%mLayerVolFracIce)%dat(1:nSnow), & ! intent(inout): volumetric fraction of ice after itertations (-)
+ !********** SUNDIALS ADDITION *****************
+ if (sundials) then
+ call computSnowDepth(&
+ dt_sub, & ! intent(in)
+ nSnow, & ! intent(in)
+ scalarSnowSublimation, & ! intent(in)
+ mLayerVolFracLiq, & ! intent(inout)
+ mLayerVolFracIce, & ! intent(inout)
+ prog_data%var(iLookPROG%mLayerTemp)%dat, & ! intent(in)
+ diag_data%var(iLookDIAG%mLayerMeltFreeze)%dat,& ! intent(in)
+ mpar_data, & ! intent(in)
+ ! output
+ tooMuchSublim, & ! intent(out): flag to denote that there was too much sublimation in a given time step
+ mLayerDepth, & ! intent(inout)
+ ! error control
+ err,message) ! intent(out): error control
+ if(err/=0)then
+ err=55
+ print*, "computSnowDepth", message
+ return
+ end if
+
+ ! process the flag for too much sublimation
+ if(tooMuchSublim)then
+ stepFailure = .true.
+ doLayerMerge = .true.
+ else
+ doLayerMerge = .false.
+ endif
+
+ ! handle special case of the step failure
+ ! NOTE: need to revert back to the previous state vector that we were happy with and reduce the time step
+ if(stepFailure)then
+ ! halve step
+ dt_sub = dtSave/2._rkind
+ ! check that the step is not tiny
+ if(dt_sub < minstep)then
+ print*,ixSolution
+ print*, 'dtSave, dt_sub', dtSave, dt_sub
+ message=trim(message)//'length of the coupled step is below the minimum step length'
+ err=20
+ return
+ endif
+
+ ! try again
+ cycle substeps
+ endif
+
+ ! end associate sublime
+
+ ! update coordinate variables
+ call calcHeight(&
+ ! input/output: data structures
+ indx_data, & ! intent(in): layer type
+ prog_data, & ! intent(inout): model variables for a local HRU
+ ! output: error control
+ err,cmessage)
+ if(err/=0)then
+ err=20
+ message=trim(message)//trim(cmessage)
+ return
+ end if
+
+ ! recompute snow depth and SWE
+ if(nSnow > 0)then
+ prog_data%var(iLookPROG%scalarSnowDepth)%dat(1) = sum( prog_data%var(iLookPROG%mLayerDepth)%dat(1:nSnow))
+ prog_data%var(iLookPROG%scalarSWE)%dat(1) = sum( (prog_data%var(iLookPROG%mLayerVolFracLiq)%dat(1:nSnow)*iden_water + &
+ prog_data%var(iLookPROG%mLayerVolFracIce)%dat(1:nSnow)*iden_ice) &
+ * prog_data%var(iLookPROG%mLayerDepth)%dat(1:nSnow) )
+ end if
+
+ ! increment fluxes
+ dt_wght = dt_sub/data_step ! define weight applied to each sub-step
+ do iVar=1,size(averageFlux_meta)
+ flux_mean%var(iVar)%dat(:) = flux_mean%var(iVar)%dat(:) + flux_data%var(averageFlux_meta(iVar)%ixParent)%dat(:)*dt_wght
+ end do
+
+ ! increment change in storage associated with the surface melt pond (kg m-2)
+ if(nSnow==0) sfcMeltPond = sfcMeltPond + prog_data%var(iLookPROG%scalarSfcMeltPond)%dat(1)
+
+ ! increment soil compression (kg m-2)
+ totalSoilCompress = totalSoilCompress + diag_data%var(iLookDIAG%scalarSoilCompress)%dat(1) ! total soil compression over whole layer (kg m-2)
+ !********** END SUNDIALS ADDITION *****************
+ else ! Non - Sundials version
+ ! * compute change in ice content of the top snow layer due to sublimation...
+ ! ---------------------------------------------------------------------------
+ ! NOTE: this is done BEFORE densification
+ if(nSnow > 0)then ! snow layers exist
+
+ ! try to remove ice from the top layer
+ iSnow=1
+
+ ! save the mass of liquid water (kg m-2)
+ massLiquid = mLayerDepth(iSnow)*mLayerVolFracLiq(iSnow)*iden_water
+
+ ! add/remove the depth of snow gained/lost by frost/sublimation (m)
+ ! NOTE: assume constant density
+ mLayerDepth(iSnow) = mLayerDepth(iSnow) + dt_sub*scalarSnowSublimation/(mLayerVolFracIce(iSnow)*iden_ice)
+
+ ! check that we did not remove the entire layer
+ if(mLayerDepth(iSnow) < verySmall)then
+ stepFailure = .true.
+ doLayerMerge = .true.
+ dt_sub = max(dtSave/2._dp, minstep)
+ cycle substeps
+ else
+ stepFailure = .false.
+ doLayerMerge = .false.
+ endif
+
+ ! update the volumetric fraction of liquid water
+ mLayerVolFracLiq(iSnow) = massLiquid / (mLayerDepth(iSnow)*iden_water)
+
+ ! no snow
+ else
+
+ ! no snow: check that sublimation is zero
+ if(abs(scalarSnowSublimation) > verySmall)then
+ message=trim(message)//'sublimation of snow has been computed when no snow exists'
+ print*, message
+ err=20; return
+ end if
+
+ end if ! (if snow layers exist)
+
+ ! end associate sublime
+
+ ! *** account for compaction and cavitation in the snowpack...
+ ! ------------------------------------------------------------
+ if(nSnow>0)then
+ call snwDensify(&
+ ! intent(in): variables
+ dt_sub, & ! intent(in): time step (s)
+ indx_data%var(iLookINDEX%nSnow)%dat(1), & ! intent(in): number of snow layers
+ prog_data%var(iLookPROG%mLayerTemp)%dat(1:nSnow), & ! intent(in): temperature of each layer (K)
+ diag_data%var(iLookDIAG%mLayerMeltFreeze)%dat(1:nSnow), & ! intent(in): volumetric melt in each layer (kg m-3)
+ ! intent(in): parameters
+ mpar_data%var(iLookPARAM%densScalGrowth)%dat(1), & ! intent(in): density scaling factor for grain growth (kg-1 m3)
+ mpar_data%var(iLookPARAM%tempScalGrowth)%dat(1), & ! intent(in): temperature scaling factor for grain growth (K-1)
+ mpar_data%var(iLookPARAM%grainGrowthRate)%dat(1), & ! intent(in): rate of grain growth (s-1)
+ mpar_data%var(iLookPARAM%densScalOvrbdn)%dat(1), & ! intent(in): density scaling factor for overburden pressure (kg-1 m3)
+ mpar_data%var(iLookPARAM%tempScalOvrbdn)%dat(1), & ! intent(in): temperature scaling factor for overburden pressure (K-1)
+ mpar_data%var(iLookPARAM%baseViscosity)%dat(1), & ! intent(in): viscosity coefficient at T=T_frz and snow density=0 (kg m-2 s)
+ ! intent(inout): state variables
+ prog_data%var(iLookPROG%mLayerDepth)%dat(1:nSnow), & ! intent(inout): depth of each layer (m)
+ prog_data%var(iLookPROG%mLayerVolFracLiq)%dat(1:nSnow), & ! intent(inout): volumetric fraction of liquid water after itertations (-)
+ prog_data%var(iLookPROG%mLayerVolFracIce)%dat(1:nSnow), & ! intent(inout): volumetric fraction of ice after itertations (-)
+ ! output: error control
+ err,cmessage) ! intent(out): error control
+ if(err/=0)then
+ err=55
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ end if
+
+ end if ! if snow layers exist
+
+ ! update coordinate variables
+ call calcHeight(&
+ ! input/output: data structures
+ indx_data, & ! intent(in): layer type
+ prog_data, & ! intent(inout): model variables for a local HRU
! output: error control
- err,cmessage) ! intent(out): error control
- if(err/=0)then; err=55; message=trim(message)//trim(cmessage); return; end if
- end if ! if snow layers exist
-
- ! update coordinate variables
- call calcHeight(&
- ! input/output: data structures
- indx_data, & ! intent(in): layer type
- prog_data, & ! intent(inout): model variables for a local HRU
- ! output: error control
- err,cmessage)
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if
-
- ! recompute snow depth and SWE
- if(nSnow > 0)then
- prog_data%var(iLookPROG%scalarSnowDepth)%dat(1) = sum( prog_data%var(iLookPROG%mLayerDepth)%dat(1:nSnow))
- prog_data%var(iLookPROG%scalarSWE)%dat(1) = sum( (prog_data%var(iLookPROG%mLayerVolFracLiq)%dat(1:nSnow)*iden_water + &
- prog_data%var(iLookPROG%mLayerVolFracIce)%dat(1:nSnow)*iden_ice) &
- * prog_data%var(iLookPROG%mLayerDepth)%dat(1:nSnow) )
- end if
-
- ! increment fluxes
- dt_wght = dt_sub/data_step ! define weight applied to each sub-step
- do iVar=1,size(averageFlux_meta)
- flux_mean%var(iVar)%dat(:) = flux_mean%var(iVar)%dat(:) + flux_data%var(averageFlux_meta(iVar)%ixParent)%dat(:)*dt_wght
- end do
+ err,cmessage)
+ if(err/=0)then
+ err=20
+ message=trim(message)//trim(cmessage)
+ print*, message
+ return
+ end if
+
+ ! recompute snow depth and SWE
+ if(nSnow > 0)then
+ prog_data%var(iLookPROG%scalarSnowDepth)%dat(1) = sum( prog_data%var(iLookPROG%mLayerDepth)%dat(1:nSnow))
+ prog_data%var(iLookPROG%scalarSWE)%dat(1) = sum( (prog_data%var(iLookPROG%mLayerVolFracLiq)%dat(1:nSnow)*iden_water + &
+ prog_data%var(iLookPROG%mLayerVolFracIce)%dat(1:nSnow)*iden_ice) &
+ * prog_data%var(iLookPROG%mLayerDepth)%dat(1:nSnow) )
+ end if
+
+ ! increment fluxes
+ dt_wght = dt_sub/data_step ! define weight applied to each sub-step
+ do iVar=1,size(averageFlux_meta)
+ flux_mean%var(iVar)%dat(:) = flux_mean%var(iVar)%dat(:) + flux_data%var(averageFlux_meta(iVar)%ixParent)%dat(:)*dt_wght
+ end do
+
+ ! increment change in storage associated with the surface melt pond (kg m-2)
+ if(nSnow==0) sfcMeltPond = sfcMeltPond + prog_data%var(iLookPROG%scalarSfcMeltPond)%dat(1)
+
+ ! increment soil compression (kg m-2)
+ totalSoilCompress = totalSoilCompress + diag_data%var(iLookDIAG%scalarSoilCompress)%dat(1) ! total soil compression over whole layer (kg m-2)
- ! increment change in storage associated with the surface melt pond (kg m-2)
- if(nSnow==0) sfcMeltPond = sfcMeltPond + prog_data%var(iLookPROG%scalarSfcMeltPond)%dat(1)
+ end if ! sundials
+ end associate sublime
- ! increment soil compression (kg m-2)
- totalSoilCompress = totalSoilCompress + diag_data%var(iLookDIAG%scalarSoilCompress)%dat(1) ! total soil compression over whole layer (kg m-2)
! ****************************************************************************************************
! *** END MAIN SOLVER ********************************************************************************
diff --git a/build/source/engine/nr_utility.f90 b/build/source/engine/nr_utility.f90
index 238f906..adc4d8b 100755
--- a/build/source/engine/nr_utility.f90
+++ b/build/source/engine/nr_utility.f90
@@ -12,23 +12,6 @@ public::arth
public::indexx
contains
- ! *************************************************************************************************
- ! * the arth function, used to build a vector of regularly spaced numbers
- ! *************************************************************************************************
-! FUNCTION arth_r(first,increment,n)
-! implicit none
-! REAL(SP), INTENT(IN) :: first,increment
-! INTEGER(I4B), INTENT(IN) :: n
-! REAL(SP), DIMENSION(n) :: arth_r
-! INTEGER(I4B) :: k
-! arth_r(1)=first
-! if(n>1)then
-! do k=2,n
-! arth_r(k) = arth_r(k-1) + increment
-! end do
-! end if
-! END FUNCTION arth_r
- ! ------------------------------------------------------------------------------------------------
FUNCTION arth_d(first,increment,n)
implicit none
REAL(rkind), INTENT(IN) :: first,increment
--
GitLab