diff --git a/build/cmake/CMakeLists.txt b/build/cmake/CMakeLists.txt
index 46a28fef72831326aea471121cf936142c3d31e8..19b411a40c222bb2246b412bf78f7d8848e4639a 100644
--- a/build/cmake/CMakeLists.txt
+++ b/build/cmake/CMakeLists.txt
@@ -2,299 +2,624 @@ cmake_minimum_required(VERSION 3.10 FATAL_ERROR)
set(PARENT_DIR ../../)
set(EXEC_DIR ${PARENT_DIR}../bin)
-
set(EXEC_NAME summa_actors)
-
project(summa_actors LANGUAGES CXX Fortran)
enable_language(C)
SET (CMAKE_Fortran_COMPILER gfortran)
include(FortranCInterface)
FortranCInterface_VERIFY(CXX)
-set(ACTORS_DIR ${PARENT_DIR}/build/source/actors)
-set(DRIVER_DIR ${PARENT_DIR}/build/source/driver)
-set(DSHARE_DIR ${PARENT_DIR}/build/source/dshare)
-set(ENGINE_DIR ${PARENT_DIR}/build/source/engine)
-set(HOOKUP_DIR ${PARENT_DIR}/build/source/hookup)
-set(NETCDF_DIR ${PARENT_DIR}/build/source/netcdf)
-set(NOAHMP_DIR ${PARENT_DIR}/build/source/noah-mp)
-set(FILE_ACCESS_DIR ${ACTORS_DIR}/file_access_actor)
-set(JOB_ACTOR_DIR ${ACTORS_DIR}/job_actor)
-set(HRU_ACTOR_DIR ${ACTORS_DIR}/hru_actor)
-set(GRU_ACTOR_DIR ${ACTORS_DIR}/gru_actor)
-
-set(SUMMA_DSHARE_DIR ${PARENT_DIR}/build/summa/build/source/dshare)
-set(SUMMA_ENGINE_DIR ${PARENT_DIR}/build/summa/build/source/engine)
-
-set(NRUTIL
- ${ENGINE_DIR}/nrtype.f90
- ${ENGINE_DIR}/f2008funcs.f90
- ${ENGINE_DIR}/nr_utility.f90)
-
-set(NRPROC
- ${ENGINE_DIR}/expIntegral.f90
- ${ENGINE_DIR}/spline_int.f90)
-
-SET(HOOKUP
- ${HOOKUP_DIR}/ascii_util.f90
- ${HOOKUP_DIR}/summaActors_FileManager.f90)
-
-SET(DATAMS
- ${DSHARE_DIR}/multiconst.f90
- ${DSHARE_DIR}/var_lookup.f90
- ${DSHARE_DIR}/data_types.f90
- ${DSHARE_DIR}/globalData.f90
- ${DSHARE_DIR}/flxMapping.f90)
-
-SET(DEPENDS_ON_FILEMANAGER
- ${SUMMA_DSHARE_DIR}/get_ixname.f90
- ${SUMMA_DSHARE_DIR}/popMetadat.f90
- ${DSHARE_DIR}/outpt_stat.f90)
-
-SET(UTILMS
- ${ENGINE_DIR}/time_utils.f90
- ${ENGINE_DIR}/mDecisions.f90
- ${ENGINE_DIR}/snow_utils.f90
- ${ENGINE_DIR}/soil_utils.f90
- ${ENGINE_DIR}/updatState.f90
- ${ENGINE_DIR}/matrixOper.f90)
-
-set(SOLVER
- ${ENGINE_DIR}/vegPhenlgy.f90
- ${ENGINE_DIR}/sundials/computSnowDepth.f90
- ${SUMMA_ENGINE_DIR}/diagn_evar.f90
- ${SUMMA_ENGINE_DIR}/stomResist.f90
- ${SUMMA_ENGINE_DIR}/groundwatr.f90
- ${SUMMA_ENGINE_DIR}/vegSWavRad.f90
- ${SUMMA_ENGINE_DIR}/vegNrgFlux.f90
- ${SUMMA_ENGINE_DIR}/ssdNrgFlux.f90
- ${SUMMA_ENGINE_DIR}/vegLiqFlux.f90
- ${SUMMA_ENGINE_DIR}/snowLiqFlx.f90
- ${SUMMA_ENGINE_DIR}/soilLiqFlx.f90
- ${SUMMA_ENGINE_DIR}/bigAquifer.f90
- ${SUMMA_ENGINE_DIR}/computFlux.f90
- ${SUMMA_ENGINE_DIR}/computResid.f90
- ${SUMMA_ENGINE_DIR}/computJacob.f90
- ${SUMMA_ENGINE_DIR}/eval8summa.f90
- ${SUMMA_ENGINE_DIR}/summaSolve.f90
- ${SUMMA_ENGINE_DIR}/systemSolv.f90
- ${SUMMA_ENGINE_DIR}/varSubstep.f90
- ${SUMMA_ENGINE_DIR}/opSplittin.f90
- ${ENGINE_DIR}/coupled_em.f90)
-
-set(INTERFACE
- ${ACTORS_DIR}/global/cppwrap_datatypes.f90
- ${ACTORS_DIR}/global/cppwrap_auxiliary.f90
- ${ACTORS_DIR}/global/cppwrap_metadata.f90)
-
-set(FILE_ACCESS_INTERFACE
- ${FILE_ACCESS_DIR}/fortran_code/output_structure.f90
- ${FILE_ACCESS_DIR}/fortran_code/cppwrap_fileAccess.f90
- ${FILE_ACCESS_DIR}/fortran_code/read_attribute.f90
- ${FILE_ACCESS_DIR}/fortran_code/read_forcing.f90
- ${FILE_ACCESS_DIR}/fortran_code/read_param.f90
- ${FILE_ACCESS_DIR}/fortran_code/read_initcond.f90
- ${FILE_ACCESS_DIR}/fortran_code/writeOutputFromOutputStructure.f90
- ${FILE_ACCESS_DIR}/fortran_code/write_to_netcdf.f90)
-
-set(JOB_INTERFACE
- ${JOB_ACTOR_DIR}/job_actor.f90)
-
-set(HRU_INTERFACE
- ${HRU_ACTOR_DIR}/fortran_code/cppwrap_hru.f90
- ${HRU_ACTOR_DIR}/fortran_code/hru_actor.f90
- ${HRU_ACTOR_DIR}/fortran_code/init_hru_actor.f90
- ${HRU_ACTOR_DIR}/fortran_code/outputStrucWrite.f90
- ${HRU_ACTOR_DIR}/fortran_code/hru_writeOutput.f90)
-
-set(GRU_INTERFACE
- ${GRU_ACTOR_DIR}/gru_actor.f90)
-
-set(PRELIM
- ${ENGINE_DIR}/conv_funcs.f90
- ${ENGINE_DIR}/sunGeomtry.f90
- ${ENGINE_DIR}/convE2Temp.f90
- ${ENGINE_DIR}/allocspaceActors.f90
- ${ENGINE_DIR}/alloc_fileAccess.f90
- ${ENGINE_DIR}/checkStruc.f90
- ${ENGINE_DIR}/childStruc.f90
- ${ENGINE_DIR}/ffile_info.f90
- ${ENGINE_DIR}/read_dimension.f90
- ${ENGINE_DIR}/read_pinit.f90
- ${ENGINE_DIR}/pOverwrite.f90
- ${ENGINE_DIR}/paramCheck.f90
- ${ENGINE_DIR}/check_icondActors.f90)
-
-set(NOAHMP
- ${NOAHMP_DIR}/module_model_constants.F
- ${NOAHMP_DIR}/module_sf_noahutl.F
- ${NOAHMP_DIR}/module_sf_noahlsm.F
- ${NOAHMP_DIR}/module_sf_noahmplsm.F)
-
-set(MODRUN
- ${SUMMA_ENGINE_DIR}/indexState.f90
- ${SUMMA_ENGINE_DIR}/getVectorz.f90
- ${SUMMA_ENGINE_DIR}/updateVars.f90
- ${SUMMA_ENGINE_DIR}/var_derive.f90
- ${ENGINE_DIR}/derivforce.f90
- ${ENGINE_DIR}/sundials/t2enthalpy.f90
- ${SUMMA_ENGINE_DIR}/snowAlbedo.f90
- ${SUMMA_ENGINE_DIR}/canopySnow.f90
- ${SUMMA_ENGINE_DIR}/tempAdjust.f90
- ${SUMMA_ENGINE_DIR}/snwCompact.f90
- ${SUMMA_ENGINE_DIR}/layerMerge.f90
- ${SUMMA_ENGINE_DIR}/layerDivide.f90
- ${SUMMA_ENGINE_DIR}/volicePack.f90
- ${SUMMA_ENGINE_DIR}/qTimeDelay.f90)
-
-set(NETCDF
- ${NETCDF_DIR}/netcdf_util.f90
- ${NETCDF_DIR}/def_output.f90
- ${NETCDF_DIR}/writeOutput.f90
- ${NETCDF_DIR}/read_icondActors.f90)
-
-set(DRIVER
- ${DRIVER_DIR}/summaActors_type.f90
- ${DRIVER_DIR}/summaActors_util.f90
- ${DRIVER_DIR}/summaActors_globalData.f90
- ${DRIVER_DIR}/SummaActors_setup.f90
- ${DRIVER_DIR}/summaActors_restart.f90
- ${DRIVER_DIR}/SummaActors_modelRun.f90
- ${DRIVER_DIR}/summaActors_alarms.f90)
-
-set(COMM_ALL
- ${NRPROC}
- ${DATAMS}
- ${INTERFACE}
- ${HOOKUP}
- ${DEPENDS_ON_FILEMANAGER}
- ${UTILMS})
-
-set(SUMMA_ALL
- ${NETCDF}
- ${PRELIM}
- ${MODRUN}
- ${SOLVER}
- ${DRIVER}
- ${JOB_INTERFACE}
- ${FILE_ACCESS_INTERFACE}
- ${HRU_INTERFACE}
- ${GRU_INTERFACE})
-
-set(ACTORS_GLOBAL
- ${ACTORS_DIR}/global/global.cpp
- ${ACTORS_DIR}/global/timing_info.cpp
- ${ACTORS_DIR}/global/settings_functions.cpp
- ${ACTORS_DIR}/global/auxiliary.cpp)
-
-set(SUMMA_ACTOR
- ${ACTORS_DIR}/summa_actor/summa_actor.cpp
- ${ACTORS_DIR}/summa_actor/summa_client.cpp
- ${ACTORS_DIR}/summa_actor/summa_server.cpp
- ${ACTORS_DIR}/summa_actor/summa_backup_server.cpp
- ${ACTORS_DIR}/summa_actor/batch/batch.cpp
- ${ACTORS_DIR}/summa_actor/batch/batch_container.cpp
- ${ACTORS_DIR}/summa_actor/client/client.cpp
- ${ACTORS_DIR}/summa_actor/client/client_container.cpp)
-
-set(HRU_ACTOR
- ${ACTORS_DIR}/hru_actor/cpp_code/hru_actor.cpp)
-
-set(GRU_ACTOR
- ${ACTORS_DIR}/gru_actor/gru_actor.cpp)
-
-set(FILE_ACCESS_ACTOR
- ${ACTORS_DIR}/file_access_actor/cpp_code/file_access_actor.cpp
- ${ACTORS_DIR}/file_access_actor/cpp_code/forcing_file_info.cpp
- ${ACTORS_DIR}/file_access_actor/cpp_code/output_container.cpp)
-
-set(JOB_ACTOR
- ${ACTORS_DIR}/job_actor/job_actor.cpp
- ${ACTORS_DIR}/job_actor/GRUinfo.cpp)
-
-set(MAIN
- ${ACTORS_DIR}/main.cpp)
-
-
+option(SUNDIALS "Use SUNDIALS" OFF)
+
+if (SUNDIALS)
+ set(ACTORS_DIR ${PARENT_DIR}/build/source/actors)
+ set(DRIVER_DIR ${PARENT_DIR}/build/source/driver)
+ set(DSHARE_DIR ${PARENT_DIR}/build/source/dshare)
+ set(ENGINE_DIR ${PARENT_DIR}/build/source/engine)
+ set(HOOKUP_DIR ${PARENT_DIR}/build/source/hookup)
+ set(NETCDF_DIR ${PARENT_DIR}/build/source/netcdf)
+ set(NOAHMP_DIR ${PARENT_DIR}/build/source/noah-mp)
+ set(FILE_ACCESS_DIR ${ACTORS_DIR}/file_access_actor)
+ set(JOB_ACTOR_DIR ${ACTORS_DIR}/job_actor)
+ set(HRU_ACTOR_DIR ${ACTORS_DIR}/hru_actor)
+ set(GRU_ACTOR_DIR ${ACTORS_DIR}/gru_actor)
+ set(SUMMA_DSHARE_DIR ${PARENT_DIR}/build/summa-sundials/build/source/dshare)
+ set(SUMMA_ENGINE_DIR ${PARENT_DIR}/build/summa-sundials/build/source/engine)
+
+ set(DIR_SUNDIALS /usr/local/sundials)
+
+ set(NRUTIL
+ ${SUMMA_ENGINE_DIR}/nrtype.f90
+ ${SUMMA_ENGINE_DIR}/f2008funcs.f90
+ ${SUMMA_ENGINE_DIR}/nr_utility.f90)
+
+ set(NRPROC
+ ${ENGINE_DIR}/expIntegral.f90
+ ${SUMMA_ENGINE_DIR}/spline_int.f90)
+
+ SET(HOOKUP
+ ${HOOKUP_DIR}/ascii_util.f90
+ ${HOOKUP_DIR}/summaActors_FileManager.f90)
+
+ SET(DATAMS
+ ${DSHARE_DIR}/multiconst.f90
+ ${SUMMA_DSHARE_DIR}/var_lookup.f90
+ ${DSHARE_DIR}/data_types.f90
+ ${DSHARE_DIR}/globalData.f90
+ ${DSHARE_DIR}/flxMapping.f90)
+
+ SET(DEPENDS_ON_FILEMANAGER
+ ${SUMMA_DSHARE_DIR}/get_ixname.f90
+ ${SUMMA_DSHARE_DIR}/popMetadat.f90
+ ${DSHARE_DIR}/outpt_stat.f90)
+
+ SET(UTILMS
+ ${ENGINE_DIR}/time_utils.f90
+ ${ENGINE_DIR}/sundials/mDecisions.f90
+ ${SUMMA_ENGINE_DIR}/snow_utils.f90
+ ${SUMMA_ENGINE_DIR}/soil_utils.f90
+ ${SUMMA_ENGINE_DIR}/soil_utilsAddSundials.f90
+ ${SUMMA_ENGINE_DIR}/updatState.f90
+ ${SUMMA_ENGINE_DIR}/updatStateSundials.f90
+ ${SUMMA_ENGINE_DIR}/matrixOper.f90)
+
+ set(SOLVER
+ ${ENGINE_DIR}/vegPhenlgy.f90
+ ${SUMMA_ENGINE_DIR}/diagn_evar.f90
+ ${SUMMA_ENGINE_DIR}/stomResist.f90
+ ${SUMMA_ENGINE_DIR}/groundwatr.f90
+ ${SUMMA_ENGINE_DIR}/vegSWavRad.f90
+ ${SUMMA_ENGINE_DIR}/vegNrgFlux.f90
+ ${SUMMA_ENGINE_DIR}/ssdNrgFlux.f90
+ ${SUMMA_ENGINE_DIR}/vegLiqFlux.f90
+ ${SUMMA_ENGINE_DIR}/snowLiqFlx.f90
+ ${SUMMA_ENGINE_DIR}/soilLiqFlx.f90
+ ${SUMMA_ENGINE_DIR}/bigAquifer.f90
+ ${SUMMA_ENGINE_DIR}/computFlux.f90
+ ${SUMMA_ENGINE_DIR}/type4IDA.f90
+ ${SUMMA_ENGINE_DIR}/tol4IDA.f90
+ ${SUMMA_ENGINE_DIR}/computEnthalpy.f90
+ ${SUMMA_ENGINE_DIR}/computHeatCap.f90
+ ${SUMMA_ENGINE_DIR}/computThermConduct.f90
+ ${SUMMA_ENGINE_DIR}/computResid.f90
+ ${SUMMA_ENGINE_DIR}/computJacob.f90
+ ${SUMMA_ENGINE_DIR}/eval8summa.f90
+ ${SUMMA_ENGINE_DIR}/summaSolve.f90
+ ${SUMMA_ENGINE_DIR}/systemSolv.f90
+ ${SUMMA_ENGINE_DIR}/computResidSundials.f90
+ ${SUMMA_ENGINE_DIR}/eval8summaSundials.f90
+ ${SUMMA_ENGINE_DIR}/computJacobSundials.f90
+ ${SUMMA_ENGINE_DIR}/computSnowDepth.f90
+ ${SUMMA_ENGINE_DIR}/summaSolveSundialsIDA.f90
+ ${SUMMA_ENGINE_DIR}/systemSolvSundials.f90
+ ${SUMMA_ENGINE_DIR}/varSubstep.f90
+ ${SUMMA_ENGINE_DIR}/varSubstepSundials.f90
+ ${SUMMA_ENGINE_DIR}/varSubstep.f90
+ ${SUMMA_ENGINE_DIR}/opSplittin.f90
+ ${ENGINE_DIR}/sundials/coupled_em.f90)
+
+ set(INTERFACE
+ ${ACTORS_DIR}/global/cppwrap_datatypes.f90
+ ${ACTORS_DIR}/global/cppwrap_auxiliary.f90
+ ${ACTORS_DIR}/global/cppwrap_metadata.f90)
+
+ set(FILE_ACCESS_INTERFACE
+ ${FILE_ACCESS_DIR}/fortran_code/output_structure.f90
+ ${FILE_ACCESS_DIR}/fortran_code/cppwrap_fileAccess.f90
+ ${FILE_ACCESS_DIR}/fortran_code/read_attribute.f90
+ ${FILE_ACCESS_DIR}/fortran_code/read_forcing.f90
+ ${FILE_ACCESS_DIR}/fortran_code/read_param.f90
+ ${FILE_ACCESS_DIR}/fortran_code/read_initcond.f90
+ ${FILE_ACCESS_DIR}/fortran_code/writeOutputFromOutputStructure.f90
+ ${FILE_ACCESS_DIR}/fortran_code/write_to_netcdf.f90)
+
+ set(JOB_INTERFACE
+ ${JOB_ACTOR_DIR}/job_actor.f90)
+
+ set(HRU_INTERFACE
+ ${HRU_ACTOR_DIR}/fortran_code/cppwrap_hru.f90
+ ${HRU_ACTOR_DIR}/fortran_code/hru_actor.f90
+ ${HRU_ACTOR_DIR}/fortran_code/init_hru_actor.f90
+ ${HRU_ACTOR_DIR}/fortran_code/outputStrucWrite.f90
+ ${HRU_ACTOR_DIR}/fortran_code/hru_writeOutput.f90)
+
+ # set(GRU_INTERFACE
+ # ${GRU_ACTOR_DIR}/gru_actor.f90)
+
+ set(PRELIM
+ ${ENGINE_DIR}/conv_funcs.f90
+ ${ENGINE_DIR}/sunGeomtry.f90
+ ${SUMMA_ENGINE_DIR}/convE2Temp.f90
+ ${ENGINE_DIR}/allocspaceActors.f90
+ ${ENGINE_DIR}/alloc_fileAccess.f90
+ ${ENGINE_DIR}/checkStruc.f90
+ ${ENGINE_DIR}/childStruc.f90
+ ${ENGINE_DIR}/ffile_info.f90
+ ${ENGINE_DIR}/read_dimension.f90
+ ${ENGINE_DIR}/read_pinit.f90
+ ${ENGINE_DIR}/pOverwrite.f90
+ ${ENGINE_DIR}/paramCheck.f90
+ ${ENGINE_DIR}/check_icondActors.f90)
+
+ set(NOAHMP
+ ${NOAHMP_DIR}/module_model_constants.F
+ ${NOAHMP_DIR}/module_sf_noahutl.F
+ ${NOAHMP_DIR}/module_sf_noahlsm.F
+ ${NOAHMP_DIR}/module_sf_noahmplsm.F)
+
+ set(MODRUN
+ ${SUMMA_ENGINE_DIR}/indexState.f90
+ ${SUMMA_ENGINE_DIR}/getVectorz.f90
+ ${SUMMA_ENGINE_DIR}/t2enthalpy.f90
+ ${SUMMA_ENGINE_DIR}/updateVars.f90
+ ${SUMMA_ENGINE_DIR}/updateVarsSundials.f90
+ ${SUMMA_ENGINE_DIR}/var_derive.f90
+ ${ENGINE_DIR}/derivforce.f90
+ ${SUMMA_ENGINE_DIR}/snowAlbedo.f90
+ ${SUMMA_ENGINE_DIR}/canopySnow.f90
+ ${SUMMA_ENGINE_DIR}/tempAdjust.f90
+ ${SUMMA_ENGINE_DIR}/snwCompact.f90
+ ${SUMMA_ENGINE_DIR}/layerMerge.f90
+ ${SUMMA_ENGINE_DIR}/layerDivide.f90
+ ${SUMMA_ENGINE_DIR}/volicePack.f90
+ ${SUMMA_ENGINE_DIR}/qTimeDelay.f90)
+
+ set(NETCDF
+ ${NETCDF_DIR}/netcdf_util.f90
+ ${NETCDF_DIR}/def_output.f90
+ ${NETCDF_DIR}/writeOutput.f90
+ ${NETCDF_DIR}/read_icondActors.f90)
+
+ set(DRIVER
+ ${DRIVER_DIR}/summaActors_type.f90
+ ${DRIVER_DIR}/summaActors_util.f90
+ ${DRIVER_DIR}/summaActors_globalData.f90
+ ${DRIVER_DIR}/SummaActors_setup.f90
+ ${DRIVER_DIR}/summaActors_restart.f90
+ ${DRIVER_DIR}/SummaActors_modelRun.f90
+ ${DRIVER_DIR}/summaActors_alarms.f90)
+
+ set(COMM_ALL
+ ${NRPROC}
+ ${DATAMS}
+ ${INTERFACE}
+ ${HOOKUP}
+ ${DEPENDS_ON_FILEMANAGER}
+ ${UTILMS})
+
+ set(SUMMA_ALL
+ ${NETCDF}
+ ${PRELIM}
+ ${MODRUN}
+ ${SOLVER}
+ ${DRIVER}
+ ${JOB_INTERFACE}
+ ${FILE_ACCESS_INTERFACE}
+ ${HRU_INTERFACE}
+ ${GRU_INTERFACE})
+
+ set(ACTORS_GLOBAL
+ ${ACTORS_DIR}/global/global.cpp
+ ${ACTORS_DIR}/global/timing_info.cpp
+ ${ACTORS_DIR}/global/message_atoms.cpp
+ ${ACTORS_DIR}/global/settings_functions.cpp
+ ${ACTORS_DIR}/global/auxiliary.cpp)
+
+ set(SUMMA_ACTOR
+ ${ACTORS_DIR}/summa_actor/summa_actor.cpp
+ ${ACTORS_DIR}/summa_actor/summa_client.cpp
+ ${ACTORS_DIR}/summa_actor/summa_server.cpp
+ ${ACTORS_DIR}/summa_actor/summa_backup_server.cpp
+ ${ACTORS_DIR}/summa_actor/batch/batch.cpp
+ ${ACTORS_DIR}/summa_actor/batch/batch_container.cpp
+ ${ACTORS_DIR}/summa_actor/client/client.cpp
+ ${ACTORS_DIR}/summa_actor/client/client_container.cpp)
+
+ set(HRU_ACTOR
+ ${ACTORS_DIR}/hru_actor/cpp_code/hru_actor.cpp)
+
+ # set(GRU_ACTOR
+ # ${ACTORS_DIR}/gru_actor/gru_actor.cpp)
+
+ set(FILE_ACCESS_ACTOR
+ ${ACTORS_DIR}/file_access_actor/cpp_code/file_access_actor.cpp
+ ${ACTORS_DIR}/file_access_actor/cpp_code/forcing_file_info.cpp
+ ${ACTORS_DIR}/file_access_actor/cpp_code/output_container.cpp)
+
+ set(JOB_ACTOR
+ ${ACTORS_DIR}/job_actor/job_actor.cpp
+ ${ACTORS_DIR}/job_actor/GRUinfo.cpp)
+
+ set(MAIN
+ ${ACTORS_DIR}/main.cpp)
+
find_package(CAF REQUIRED)
find_package(netCDF REQUIRED)
find_package(LAPACK REQUIRED)
-
-add_library(SUMMA_NOAHMP OBJECT
- ${NOAHMP}
- ${NRUTIL})
- target_compile_options(SUMMA_NOAHMP PRIVATE
- -g -O3 -ffree-form -ffree-line-length-none -fmax-errors=0 -fPIC -Wfatal-errors)
-
-add_library(SUMMA_COMM OBJECT
- ${COMM_ALL})
- target_compile_options(SUMMA_COMM PRIVATE
- -g -O3 -ffree-line-length-none -fmax-errors=0 -fPIC -Wfatal-errors)
- target_include_directories(SUMMA_COMM PRIVATE
- "/usr/include"
- "$ENV{EBROOTNETCDFMINFORTRAN}/include"
- ${netCDF_INCLUDES}
- ${LAPACK_INCLUDES})
- target_link_libraries(SUMMA_COMM PUBLIC
- "$ENV{EBROOTNETCDFMINFORTRAN}/lib64"
+ find_package(SUNDIALS REQUIRED)
+ link_directories(${DIR_SUNDIALS}/lib)
+
+ add_library(SUMMA_NOAHMP OBJECT
+ ${NOAHMP}
+ ${NRUTIL})
+ target_compile_options(SUMMA_NOAHMP PRIVATE
+ -g -O3 -ffree-form -ffree-line-length-none -fmax-errors=0 -fPIC -Wfatal-errors)
+ add_library(SUMMA_COMM OBJECT
+ ${COMM_ALL})
+ target_compile_options(SUMMA_COMM PRIVATE
+ -g -O3 -ffree-line-length-none -fmax-errors=0 -fPIC -Wfatal-errors)
+ target_include_directories(SUMMA_COMM PRIVATE
+ "/usr/include"
+ # "$ENV{EBROOTNETCDFMINFORTRAN}/include"
+ "${DIR_SUNDIALS}/include"
+ "${DIR_SUNDIALS}/fortran"
+ ${netCDF_INCLUDES}
+ ${LAPACK_INCLUDES})
+ target_link_libraries(SUMMA_COMM PUBLIC
+ # "$ENV{EBROOTNETCDFMINFORTRAN}/lib64"
+ -lsundials_fnvecmanyvector_mod
+ -lsundials_fida_mod
+ -lsundials_fnvecserial_mod
+ -lsundials_fsunlinsoldense_mod
+ -lsundials_fsunmatrixdense_mod
+ ${netCDF_LIBRARIES}
+ ${LAPACK_LIBRARIES}
+ SUMMA_NOAHMP)
+
+ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${EXEC_DIR})
+
+
+
+ add_library(summa SHARED
+ ${SUMMA_ALL})
+ target_compile_options(summa PRIVATE
+ -g -O3 -ffree-line-length-none -fmax-errors=0 -fPIC -Wfatal-errors)
+ target_include_directories(summa PUBLIC
+ "/usr/include"
+ "${DIR_SUNDIALS}/include"
+ "${DIR_SUNDIALS}/fortran"
+ # "$ENV{EBROOTNETCDFMINFORTRAN}/include"
+ ${netCDF_INCLUDES}
+ ${LAPACK_INCLUDES})
+ target_link_libraries(summa PUBLIC
+ # "$ENV{EBROOTNETCDFMINFORTRAN}/lib64"
+ -lsundials_fnvecmanyvector_mod
+ -lsundials_fida_mod
+ -lsundials_fnvecserial_mod
+ -lsundials_fsunlinsoldense_mod
+ -lsundials_fsunmatrixdense_mod
+ ${netCDF_LIBRARIES}
+ ${LAPACK_LIBRARIES}
+ SUMMA_COMM
+ SUMMA_NOAHMP
+ -lnetcdff)
+
+ set(CMAKE_CXX_FLAGS "-g -O3 -Wfatal-errors -std=c++17")
+
+
+
+ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${EXEC_DIR})
+
+ add_executable(${EXEC_NAME}
+ ${ACTORS_GLOBAL}
+ ${HRU_ACTOR}
+ ${FILE_ACCESS_ACTOR}
+ ${JOB_ACTOR}
+ ${SUMMA_ACTOR}
+ ${SUMMA_CLIENT}
+ ${SUMMA_SERVER}
+ ${MAIN})
+ set_property(TARGET ${EXEC_NAME} PROPERTY LINKER_LANGUAGE Fortran)
+ target_include_directories(${EXEC_NAME} PUBLIC
+ ${CAF_INCLUDES}
+ ${netCDF_INCLUDES}
+ ${LAPACK_INCLUDES}
+ "${DIR_SUNDIALS}/include"
+ "${PARENT_DIR}/build/includes/global"
+ "${PARENT_DIR}/build/includes/summa_actor"
+ "${PARENT_DIR}/build/includes/gru_actor"
+ "${PARENT_DIR}/build/includes/job_actor"
+ "${PARENT_DIR}/build/includes/file_access_actor"
+ "${PARENT_DIR}/build/includes/hru_actor")
+ target_link_libraries( ${EXEC_NAME}
+ ${CAF_LIBRARIES}
${netCDF_LIBRARIES}
${LAPACK_LIBRARIES}
- SUMMA_NOAHMP)
-
-set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${EXEC_DIR})
-
-
-
-add_library(summa SHARED
- ${SUMMA_ALL})
-target_compile_options(summa PRIVATE
- -g -O3 -ffree-line-length-none -fmax-errors=0 -fPIC -Wfatal-errors)
-target_include_directories(summa PUBLIC
- "/usr/include"
- "$ENV{EBROOTNETCDFMINFORTRAN}/include"
- ${netCDF_INCLUDES}
- ${LAPACK_INCLUDES})
-target_link_libraries(summa PUBLIC
- "$ENV{EBROOTNETCDFMINFORTRAN}/lib64"
+ -lcaf_core
+ -lcaf_io
+ summa
+ -lnetcdff
+ -lsundials_fnvecmanyvector_mod
+ -lsundials_fida_mod
+ -lsundials_fnvecserial_mod
+ -lsundials_fsunlinsoldense_mod
+ -lsundials_fsunmatrixdense_mod)
+
+else()
+ set(ACTORS_DIR ${PARENT_DIR}/build/source/actors)
+ set(DRIVER_DIR ${PARENT_DIR}/build/source/driver)
+ set(DSHARE_DIR ${PARENT_DIR}/build/source/dshare)
+ set(ENGINE_DIR ${PARENT_DIR}/build/source/engine)
+ set(HOOKUP_DIR ${PARENT_DIR}/build/source/hookup)
+ set(NETCDF_DIR ${PARENT_DIR}/build/source/netcdf)
+ set(NOAHMP_DIR ${PARENT_DIR}/build/source/noah-mp)
+ set(FILE_ACCESS_DIR ${ACTORS_DIR}/file_access_actor)
+ set(JOB_ACTOR_DIR ${ACTORS_DIR}/job_actor)
+ set(HRU_ACTOR_DIR ${ACTORS_DIR}/hru_actor)
+ set(GRU_ACTOR_DIR ${ACTORS_DIR}/gru_actor)
+ set(SUMMA_DSHARE_DIR ${PARENT_DIR}/build/summa/build/source/dshare)
+ set(SUMMA_ENGINE_DIR ${PARENT_DIR}/build/summa/build/source/engine)
+
+ set(NRUTIL
+ ${ENGINE_DIR}/nrtype.f90
+ ${ENGINE_DIR}/f2008funcs.f90
+ ${ENGINE_DIR}/nr_utility.f90)
+
+ set(NRPROC
+ ${ENGINE_DIR}/expIntegral.f90
+ ${ENGINE_DIR}/spline_int.f90)
+
+ SET(HOOKUP
+ ${HOOKUP_DIR}/ascii_util.f90
+ ${HOOKUP_DIR}/summaActors_FileManager.f90)
+
+ SET(DATAMS
+ ${DSHARE_DIR}/multiconst.f90
+ ${DSHARE_DIR}/var_lookup.f90
+ ${DSHARE_DIR}/data_types.f90
+ ${DSHARE_DIR}/globalData.f90
+ ${DSHARE_DIR}/flxMapping.f90)
+
+ SET(DEPENDS_ON_FILEMANAGER
+ ${SUMMA_DSHARE_DIR}/get_ixname.f90
+ ${SUMMA_DSHARE_DIR}/popMetadat.f90
+ ${DSHARE_DIR}/outpt_stat.f90)
+
+ SET(UTILMS
+ ${ENGINE_DIR}/time_utils.f90
+ ${ENGINE_DIR}/mDecisions.f90
+ ${ENGINE_DIR}/snow_utils.f90
+ ${ENGINE_DIR}/soil_utils.f90
+ ${ENGINE_DIR}/updatState.f90
+ ${ENGINE_DIR}/matrixOper.f90)
+
+ set(SOLVER
+ ${ENGINE_DIR}/vegPhenlgy.f90
+ ${ENGINE_DIR}/sundials/computSnowDepth.f90
+ ${SUMMA_ENGINE_DIR}/diagn_evar.f90
+ ${SUMMA_ENGINE_DIR}/stomResist.f90
+ ${SUMMA_ENGINE_DIR}/groundwatr.f90
+ ${SUMMA_ENGINE_DIR}/vegSWavRad.f90
+ ${SUMMA_ENGINE_DIR}/vegNrgFlux.f90
+ ${SUMMA_ENGINE_DIR}/ssdNrgFlux.f90
+ ${SUMMA_ENGINE_DIR}/vegLiqFlux.f90
+ ${SUMMA_ENGINE_DIR}/snowLiqFlx.f90
+ ${SUMMA_ENGINE_DIR}/soilLiqFlx.f90
+ ${SUMMA_ENGINE_DIR}/bigAquifer.f90
+ ${SUMMA_ENGINE_DIR}/computFlux.f90
+ ${SUMMA_ENGINE_DIR}/computResid.f90
+ ${SUMMA_ENGINE_DIR}/computJacob.f90
+ ${SUMMA_ENGINE_DIR}/eval8summa.f90
+ ${SUMMA_ENGINE_DIR}/summaSolve.f90
+ ${SUMMA_ENGINE_DIR}/systemSolv.f90
+ ${SUMMA_ENGINE_DIR}/varSubstep.f90
+ ${SUMMA_ENGINE_DIR}/opSplittin.f90
+ ${ENGINE_DIR}/coupled_em.f90)
+
+ set(INTERFACE
+ ${ACTORS_DIR}/global/cppwrap_datatypes.f90
+ ${ACTORS_DIR}/global/cppwrap_auxiliary.f90
+ ${ACTORS_DIR}/global/cppwrap_metadata.f90)
+
+ set(FILE_ACCESS_INTERFACE
+ ${FILE_ACCESS_DIR}/fortran_code/output_structure.f90
+ ${FILE_ACCESS_DIR}/fortran_code/cppwrap_fileAccess.f90
+ ${FILE_ACCESS_DIR}/fortran_code/read_attribute.f90
+ ${FILE_ACCESS_DIR}/fortran_code/read_forcing.f90
+ ${FILE_ACCESS_DIR}/fortran_code/read_param.f90
+ ${FILE_ACCESS_DIR}/fortran_code/read_initcond.f90
+ ${FILE_ACCESS_DIR}/fortran_code/writeOutputFromOutputStructure.f90
+ ${FILE_ACCESS_DIR}/fortran_code/write_to_netcdf.f90)
+
+ set(JOB_INTERFACE
+ ${JOB_ACTOR_DIR}/job_actor.f90)
+
+ set(HRU_INTERFACE
+ ${HRU_ACTOR_DIR}/fortran_code/cppwrap_hru.f90
+ ${HRU_ACTOR_DIR}/fortran_code/hru_actor.f90
+ ${HRU_ACTOR_DIR}/fortran_code/init_hru_actor.f90
+ ${HRU_ACTOR_DIR}/fortran_code/outputStrucWrite.f90
+ ${HRU_ACTOR_DIR}/fortran_code/hru_writeOutput.f90)
+
+ set(GRU_INTERFACE
+ ${GRU_ACTOR_DIR}/gru_actor.f90)
+
+ set(PRELIM
+ ${ENGINE_DIR}/conv_funcs.f90
+ ${ENGINE_DIR}/sunGeomtry.f90
+ ${ENGINE_DIR}/convE2Temp.f90
+ ${ENGINE_DIR}/allocspaceActors.f90
+ ${ENGINE_DIR}/alloc_fileAccess.f90
+ ${ENGINE_DIR}/checkStruc.f90
+ ${ENGINE_DIR}/childStruc.f90
+ ${ENGINE_DIR}/ffile_info.f90
+ ${ENGINE_DIR}/read_dimension.f90
+ ${ENGINE_DIR}/read_pinit.f90
+ ${ENGINE_DIR}/pOverwrite.f90
+ ${ENGINE_DIR}/paramCheck.f90
+ ${ENGINE_DIR}/check_icondActors.f90)
+
+ set(NOAHMP
+ ${NOAHMP_DIR}/module_model_constants.F
+ ${NOAHMP_DIR}/module_sf_noahutl.F
+ ${NOAHMP_DIR}/module_sf_noahlsm.F
+ ${NOAHMP_DIR}/module_sf_noahmplsm.F)
+
+ set(MODRUN
+ ${SUMMA_ENGINE_DIR}/indexState.f90
+ ${SUMMA_ENGINE_DIR}/getVectorz.f90
+ ${SUMMA_ENGINE_DIR}/updateVars.f90
+ ${SUMMA_ENGINE_DIR}/var_derive.f90
+ ${ENGINE_DIR}/derivforce.f90
+ ${ENGINE_DIR}/sundials/t2enthalpy.f90
+ ${SUMMA_ENGINE_DIR}/snowAlbedo.f90
+ ${SUMMA_ENGINE_DIR}/canopySnow.f90
+ ${SUMMA_ENGINE_DIR}/tempAdjust.f90
+ ${SUMMA_ENGINE_DIR}/snwCompact.f90
+ ${SUMMA_ENGINE_DIR}/layerMerge.f90
+ ${SUMMA_ENGINE_DIR}/layerDivide.f90
+ ${SUMMA_ENGINE_DIR}/volicePack.f90
+ ${SUMMA_ENGINE_DIR}/qTimeDelay.f90)
+
+ set(NETCDF
+ ${NETCDF_DIR}/netcdf_util.f90
+ ${NETCDF_DIR}/def_output.f90
+ ${NETCDF_DIR}/writeOutput.f90
+ ${NETCDF_DIR}/read_icondActors.f90)
+
+ set(DRIVER
+ ${DRIVER_DIR}/summaActors_type.f90
+ ${DRIVER_DIR}/summaActors_util.f90
+ ${DRIVER_DIR}/summaActors_globalData.f90
+ ${DRIVER_DIR}/SummaActors_setup.f90
+ ${DRIVER_DIR}/summaActors_restart.f90
+ ${DRIVER_DIR}/SummaActors_modelRun.f90
+ ${DRIVER_DIR}/summaActors_alarms.f90)
+
+ set(COMM_ALL
+ ${NRPROC}
+ ${DATAMS}
+ ${INTERFACE}
+ ${HOOKUP}
+ ${DEPENDS_ON_FILEMANAGER}
+ ${UTILMS})
+
+ set(SUMMA_ALL
+ ${NETCDF}
+ ${PRELIM}
+ ${MODRUN}
+ ${SOLVER}
+ ${DRIVER}
+ ${JOB_INTERFACE}
+ ${FILE_ACCESS_INTERFACE}
+ ${HRU_INTERFACE}
+ ${GRU_INTERFACE})
+
+ set(ACTORS_GLOBAL
+ ${ACTORS_DIR}/global/global.cpp
+ ${ACTORS_DIR}/global/timing_info.cpp
+ ${ACTORS_DIR}/global/message_atoms.cpp
+ ${ACTORS_DIR}/global/settings_functions.cpp
+ ${ACTORS_DIR}/global/auxiliary.cpp)
+
+ set(SUMMA_ACTOR
+ ${ACTORS_DIR}/summa_actor/summa_actor.cpp
+ ${ACTORS_DIR}/summa_actor/summa_client.cpp
+ ${ACTORS_DIR}/summa_actor/summa_server.cpp
+ ${ACTORS_DIR}/summa_actor/summa_backup_server.cpp
+ ${ACTORS_DIR}/summa_actor/batch/batch.cpp
+ ${ACTORS_DIR}/summa_actor/batch/batch_container.cpp
+ ${ACTORS_DIR}/summa_actor/client/client.cpp
+ ${ACTORS_DIR}/summa_actor/client/client_container.cpp)
+
+ set(HRU_ACTOR
+ ${ACTORS_DIR}/hru_actor/cpp_code/hru_actor.cpp)
+
+ set(GRU_ACTOR
+ ${ACTORS_DIR}/gru_actor/gru_actor.cpp)
+
+ set(FILE_ACCESS_ACTOR
+ ${ACTORS_DIR}/file_access_actor/cpp_code/file_access_actor.cpp
+ ${ACTORS_DIR}/file_access_actor/cpp_code/forcing_file_info.cpp
+ ${ACTORS_DIR}/file_access_actor/cpp_code/output_container.cpp)
+
+ set(JOB_ACTOR
+ ${ACTORS_DIR}/job_actor/job_actor.cpp
+ ${ACTORS_DIR}/job_actor/GRUinfo.cpp)
+
+ set(MAIN
+ ${ACTORS_DIR}/main.cpp)
+
+
+ find_package(CAF REQUIRED)
+ find_package(netCDF REQUIRED)
+ find_package(LAPACK REQUIRED)
+
+ add_library(SUMMA_NOAHMP OBJECT
+ ${NOAHMP}
+ ${NRUTIL})
+ target_compile_options(SUMMA_NOAHMP PRIVATE
+ -g -O3 -ffree-form -ffree-line-length-none -fmax-errors=0 -fPIC -Wfatal-errors)
+
+ add_library(SUMMA_COMM OBJECT
+ ${COMM_ALL})
+ target_compile_options(SUMMA_COMM PRIVATE
+ -g -O3 -ffree-line-length-none -fmax-errors=0 -fPIC -Wfatal-errors)
+ target_include_directories(SUMMA_COMM PRIVATE
+ "/usr/include"
+ # "$ENV{EBROOTNETCDFMINFORTRAN}/include"
+ ${netCDF_INCLUDES}
+ ${LAPACK_INCLUDES})
+ target_link_libraries(SUMMA_COMM PUBLIC
+ # "$ENV{EBROOTNETCDFMINFORTRAN}/lib64"
+ ${netCDF_LIBRARIES}
+ ${LAPACK_LIBRARIES}
+ SUMMA_NOAHMP)
+
+ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${EXEC_DIR})
+
+
+
+ add_library(summa SHARED
+ ${SUMMA_ALL})
+ target_compile_options(summa PRIVATE
+ -g -O3 -ffree-line-length-none -fmax-errors=0 -fPIC -Wfatal-errors)
+ target_include_directories(summa PUBLIC
+ "/usr/include"
+ # "$ENV{EBROOTNETCDFMINFORTRAN}/include"
+ ${netCDF_INCLUDES}
+ ${LAPACK_INCLUDES})
+ target_link_libraries(summa PUBLIC
+ # "$ENV{EBROOTNETCDFMINFORTRAN}/lib64"
+ ${netCDF_LIBRARIES}
+ ${LAPACK_LIBRARIES}
+ SUMMA_COMM
+ SUMMA_NOAHMP
+ -lnetcdff)
+
+ set(CMAKE_CXX_FLAGS "-g -O3 -Wfatal-errors -std=c++17")
+
+
+
+ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${EXEC_DIR})
+
+ add_executable(${EXEC_NAME}
+ ${ACTORS_GLOBAL}
+ ${HRU_ACTOR}
+ ${GRU_ACTOR}
+ ${FILE_ACCESS_ACTOR}
+ ${JOB_ACTOR}
+ ${SUMMA_ACTOR}
+ ${SUMMA_CLIENT}
+ ${SUMMA_SERVER}
+ ${MAIN})
+ set_property(TARGET ${EXEC_NAME} PROPERTY LINKER_LANGUAGE Fortran)
+ target_include_directories(${EXEC_NAME} PUBLIC
+ ${CAF_INCLUDES}
+ ${netCDF_INCLUDES}
+ ${LAPACK_INCLUDES}
+ "${PARENT_DIR}/build/includes/global"
+ "${PARENT_DIR}/build/includes/summa_actor"
+ "${PARENT_DIR}/build/includes/gru_actor"
+ "${PARENT_DIR}/build/includes/job_actor"
+ "${PARENT_DIR}/build/includes/file_access_actor"
+ "${PARENT_DIR}/build/includes/hru_actor")
+ target_link_libraries( ${EXEC_NAME}
+ ${CAF_LIBRARIES}
${netCDF_LIBRARIES}
${LAPACK_LIBRARIES}
- SUMMA_COMM
- SUMMA_NOAHMP
+ -lcaf_core
+ -lcaf_io
+ summa
-lnetcdff)
-
-set(CMAKE_CXX_FLAGS "-g -O3 -Wfatal-errors -std=c++17")
-
-
-
-set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${EXEC_DIR})
-
-add_executable(${EXEC_NAME}
- ${ACTORS_GLOBAL}
- ${HRU_ACTOR}
- ${GRU_ACTOR}
- ${FILE_ACCESS_ACTOR}
- ${JOB_ACTOR}
- ${SUMMA_ACTOR}
- ${SUMMA_CLIENT}
- ${SUMMA_SERVER}
- ${MAIN})
- set_property(TARGET ${EXEC_NAME} PROPERTY LINKER_LANGUAGE Fortran)
- target_include_directories(${EXEC_NAME} PUBLIC
- ${CAF_INCLUDES}
- ${netCDF_INCLUDES}
- ${LAPACK_INCLUDES}
- "${PARENT_DIR}/build/includes/global"
- "${PARENT_DIR}/build/includes/summa_actor"
- "${PARENT_DIR}/build/includes/gru_actor"
- "${PARENT_DIR}/build/includes/job_actor"
- "${PARENT_DIR}/build/includes/file_access_actor"
- "${PARENT_DIR}/build/includes/hru_actor")
- target_link_libraries( ${EXEC_NAME}
- ${CAF_LIBRARIES}
- ${netCDF_LIBRARIES}
- ${LAPACK_LIBRARIES}
- -lcaf_core
- -lcaf_io
- summa
- -lnetcdff
- )
\ No newline at end of file
+endif()
\ No newline at end of file
diff --git a/build/includes/global/message_atoms.hpp b/build/includes/global/message_atoms.hpp
index 35c2c8f1e5543f6a3cafab102b1dd8adc81fd837..cc7bea341e3e5bca4801e9117e499e60f6e7f22d 100644
--- a/build/includes/global/message_atoms.hpp
+++ b/build/includes/global/message_atoms.hpp
@@ -6,6 +6,20 @@
#include "client/client_container.hpp"
#include <vector>
#include "settings_functions.hpp"
+#include "caf/all.hpp"
+
+enum class hru_error : uint8_t {
+ run_physics_unhandleable = 1,
+ run_physics_infeasible_state = 2,
+};
+
+std::string to_string(hru_error err);
+bool from_string(caf::string_view in, hru_error& out);
+bool from_integer(uint8_t in, hru_error& out);
+template<class Inspector>
+bool inspect(Inspector& f, hru_error& x) {
+ return caf::default_enum_inspect(f, x);
+}
CAF_BEGIN_TYPE_ID_BLOCK(summa, first_custom_type_id)
// Sender: job_actor
@@ -187,5 +201,10 @@ CAF_BEGIN_TYPE_ID_BLOCK(summa, first_custom_type_id)
CAF_ADD_TYPE_ID(summa, (std::optional<caf::strong_actor_ptr>))
+ // error types
+ CAF_ADD_TYPE_ID(summa, (hru_error))
+
+
+CAF_END_TYPE_ID_BLOCK(summa)
-CAF_END_TYPE_ID_BLOCK(summa)
\ No newline at end of file
+CAF_ERROR_CODE_ENUM(hru_error)
\ No newline at end of file
diff --git a/build/includes/job_actor/job_actor.hpp b/build/includes/job_actor/job_actor.hpp
index 058e2b12341a268977035d9061fd5ce5f770dfa8..7c377280d808d9fc6f977e05226ce41c488456fa 100644
--- a/build/includes/job_actor/job_actor.hpp
+++ b/build/includes/job_actor/job_actor.hpp
@@ -20,7 +20,7 @@ struct job_state {
// Variables for GRU monitoring
int dt_init_start_factor = 1; // Initial Factor for dt_init (coupled_em)
- int max_run_attempts = 3; // Max number of attemtps to solve a GRU
+ int max_run_attempts = 1; // Max number of attemtps to solve a GRU
std::vector<GRUinfo*> gru_list; // List of all GRUs under this job actor
int num_gru_done = 0; // The number of GRUs that have completed
int gru_init = 0; // Number of GRUs initalized
diff --git a/build/source/actors/file_access_actor/cpp_code/file_access_actor.cpp b/build/source/actors/file_access_actor/cpp_code/file_access_actor.cpp
index fab7c20dce3f092f5c34824e3684bf9390dbb449..870a690e17d8f96ed0ce594fd8c10980462f2e9a 100644
--- a/build/source/actors/file_access_actor/cpp_code/file_access_actor.cpp
+++ b/build/source/actors/file_access_actor/cpp_code/file_access_actor.cpp
@@ -15,7 +15,7 @@ namespace caf {
behavior file_access_actor(stateful_actor<file_access_state>* self, int start_gru, int num_gru,
File_Access_Actor_Settings file_access_actor_settings, actor parent) {
aout(self) << "\n----------File_Access_Actor Started----------\n";
-
+
// Set Up timing Info we wish to track
self->state.file_access_timing = TimingInfo();
self->state.file_access_timing.addTimePoint("read_duration");
diff --git a/build/source/actors/file_access_actor/fortran_code/cppwrap_fileAccess.f90 b/build/source/actors/file_access_actor/fortran_code/cppwrap_fileAccess.f90
index 135544ab73a775552971b4fbbe51013c79d83e27..e3ccb081f2c30d831aab9c997fb393e9db40ccd7 100644
--- a/build/source/actors/file_access_actor/fortran_code/cppwrap_fileAccess.f90
+++ b/build/source/actors/file_access_actor/fortran_code/cppwrap_fileAccess.f90
@@ -10,14 +10,15 @@ module cppwrap_fileAccess
implicit none
+ public::read_pinit_C
+ public::read_vegitationTables
+ public::FileAccessActor_DeallocateStructures
public::initFailedHRUTracker
contains
-! Read in the inital parameters, from the txt files that are give to summa as input
-! LocalParamInfo.txt
-! BasinParamInfo.txt
+! Read in the inital parameters, from the txt files that are give to summa as input LocalParamInfo.txt BasinParamInfo.txt
subroutine read_pinit_C(err) bind(C, name='read_pinit_C')
USE globalData,only:localParFallback ! local column default parameters
USE globalData,only:basinParFallback ! basin-average default parameters
diff --git a/build/source/actors/global/message_atoms.cpp b/build/source/actors/global/message_atoms.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..27ecd9aa34c1a832fe500409549dfef32559a4e5
--- /dev/null
+++ b/build/source/actors/global/message_atoms.cpp
@@ -0,0 +1,44 @@
+#include "caf/all.hpp"
+#include "message_atoms.hpp"
+
+// enum class hru_error : uint8_t {
+// run_physics_unhandleable = 1,
+// run_physics_infeasible_state = 2,
+// };
+
+std::string to_string(hru_error err) {
+ switch(err) {
+ case hru_error::run_physics_unhandleable:
+ return "run_physics_unhandleable";
+ case hru_error::run_physics_infeasible_state:
+ return "run_physics_infeasible_state";
+ default:
+ return "unknown";
+ }
+}
+
+bool from_string(caf::string_view in, hru_error& out) {
+ if (in == "run_physics_unhandleable") {
+ out = hru_error::run_physics_unhandleable;
+ return true;
+ }
+ if (in == "run_physics_infeasible_state") {
+ out = hru_error::run_physics_infeasible_state;
+ return true;
+ }
+ return false;
+}
+
+bool from_integer(uint8_t in, hru_error& out) {
+ switch(in) {
+ case 1:
+ out = hru_error::run_physics_unhandleable;
+ return true;
+ case 2:
+ out = hru_error::run_physics_infeasible_state;
+ return true;
+ default:
+ return false;
+ }
+}
+
diff --git a/build/source/actors/hru_actor/cpp_code/hru_actor.cpp b/build/source/actors/hru_actor/cpp_code/hru_actor.cpp
index 28bf0580f8ab01454934250d3d1e0d8bc03583bb..1c407e7ce8f1187dc47c85ad223d0c5f0e11c296 100644
--- a/build/source/actors/hru_actor/cpp_code/hru_actor.cpp
+++ b/build/source/actors/hru_actor/cpp_code/hru_actor.cpp
@@ -139,6 +139,11 @@ behavior hru_actor(stateful_actor<hru_state>* self, int refGRU, int indxGRU,
aout(self) << "Error: HRU_Actor - Run_HRU - HRU = " << self->state.indxHRU <<
" - indxGRU = " << self->state.indxGRU << " - refGRU = "<< self->state.refGRU << std::endl;
aout(self) << "Error = " << err << "\n";
+
+ self->send(self->state.parent, run_failure_v, self, self->state.indxGRU, err);
+ // self->quit(hru_error::run_hru_unhandleable);
+ // caf::exit_reason
+ // self->down_msg(hru_error::run_physics_unhandleable);
self->quit();
}
@@ -376,6 +381,7 @@ int Run_HRU(stateful_actor<hru_state>* self) {
return 0;
}
+
void printOutput(stateful_actor<hru_state>* self) {
aout(self) << self->state.refGRU << " - Timestep = " << self->state.timestep << std::endl;
}
diff --git a/build/source/actors/job_actor/job_actor.cpp b/build/source/actors/job_actor/job_actor.cpp
index d9e3999d9b83fd716700e5a39cf13f81d172eea3..e89f16337244c260ac222dfd0b97735f5f01b717 100644
--- a/build/source/actors/job_actor/job_actor.cpp
+++ b/build/source/actors/job_actor/job_actor.cpp
@@ -18,6 +18,22 @@ behavior job_actor(stateful_actor<job_state>* self, int start_gru, int num_gru,
File_Access_Actor_Settings file_access_actor_settings, Job_Actor_Settings job_actor_settings,
HRU_Actor_Settings hru_actor_settings, caf::actor parent) {
+ self->set_down_handler([=](const down_msg& dm) {
+ aout(self) << "\n\n ********** DOWN HANDLER ********** \n";
+ aout(self) << "Lost Connection With A Connected Actor\n";
+ aout(self) << "Reason: " << to_string(dm.reason) << "\n";
+ });
+
+ self->set_error_handler([=](const error& err) {
+ aout(self) << "\n\n ********** ERROR HANDLER ********** \n";
+ aout(self) << "Error: " << to_string(err) << "\n";
+ });
+
+ self->set_exit_handler([=](const exit_msg& em) {
+ aout(self) << "\n\n ********** EXIT HANDLER ********** \n";
+ aout(self) << "Exit Reason: " << to_string(em.reason) << "\n";
+ });
+
// Timing Information
self->state.job_timing = TimingInfo();
self->state.job_timing.addTimePoint("total_duration");
@@ -118,20 +134,21 @@ behavior job_actor(stateful_actor<job_state>* self, int start_gru, int num_gru,
},
[=](run_failure, caf::actor actorRef, int indx_gru, int err) {
+
aout(self) << "GRU:" << self->state.gru_list[indx_gru - 1]->getRefGRU()
<< "indx_gru = " << indx_gru << "Failed \n"
<< "Will have to wait until all GRUs are done before it can be re-tried\n";
self->state.num_gru_failed++;
- self->state.num_gru_done++;
self->state.gru_list[indx_gru - 1]->updateFailed();
// Let the file_access_actor know this actor failed
self->send(self->state.file_access_actor, run_failure_v, indx_gru);
// check if we are the last hru to complete
- if (self->state.num_gru_done >= self->state.num_gru) {
- restartFailures(self);
+ if (self->state.num_gru_done + self->state.num_gru_failed >= self->state.num_gru) {
+ // restartFailures(self);
+ self->quit();
}
},
@@ -244,6 +261,8 @@ void restartFailures(stateful_actor<job_state>* self) {
self->state.num_gru = self->state.num_gru_failed;
self->state.num_gru_failed = 0;
self->state.num_gru_done = 0;
+
+
for(auto gru : self->state.gru_list) {
if (gru->isFailed() && !gru->isMaxAttemptsReached()) {
gru->updateFailed();
@@ -254,6 +273,13 @@ void restartFailures(stateful_actor<job_state>* self) {
gru->updateGRU(newGRU);
gru->updateCurrentAttempt();
self->send(gru->getActor(), dt_init_factor_v, gru->getDt_init());
+ } else {
+ // Max attempts reached, so we are done with this GRU
+ self->state.gru_list[gru->getIndxGRU() - 1]->doneRun(-1, -1, -1, -1, -1);
+ if (self->state.job_actor_settings.output_csv) {
+ self->state.gru_list[gru->getIndxGRU() - 1]->writeSuccess(self->state.success_output_file, self->state.hostname);
+ }
+ self->state.num_gru_done++;
}
}
}
diff --git a/build/source/dshare/globalData.f90 b/build/source/dshare/globalData.f90
index 0c92da5aaaf6b2bf1cdc022df62e9d0561338ded..41a69a6099175ecb7821c60be11e14e22f5f4fd7 100755
--- a/build/source/dshare/globalData.f90
+++ b/build/source/dshare/globalData.f90
@@ -268,6 +268,7 @@ MODULE globalData
logical(lgt),save,public :: globalPrintFlag=.false. ! flag to compute the Jacobian
integer(i4b),save,public :: chunksize=1024 ! chunk size for the netcdf read/write
integer(i4b),save,public :: outputPrecision=nf90_double ! variable type
+ integer(i4b),save,public :: outputCompressionLevel=4 ! output netcdf file deflate level: 0-9. 0 is no compression.
! define result from the time calls
integer(i4b),dimension(8),save,public :: startInit,endInit ! date/time for the start and end of the initialization
diff --git a/build/source/engine/sundials/coupled_em.f90 b/build/source/engine/sundials/coupled_em.f90
index da96f265525e79c42bcdeda31801d5d2df1f5b03..eb8022cbade0cd100cab649cb6ba7df643e7d631 100755
--- a/build/source/engine/sundials/coupled_em.f90
+++ b/build/source/engine/sundials/coupled_em.f90
@@ -203,6 +203,7 @@ subroutine coupled_em(&
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.
+ logical(lgt) :: checkMassBalance ! flag to check the mass balance
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
@@ -236,7 +237,6 @@ subroutine coupled_em(&
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=.true.
logical(lgt) :: tooMuchSublim ! flag to denote that there was too much sublimation in a given time step
@@ -914,56 +914,49 @@ subroutine coupled_em(&
mLayerVolFracIce => prog_data%var(iLookPROG%mLayerVolFracIce)%dat, & ! volumetric fraction of ice in the snow+soil domain (-)
mLayerVolFracLiq => prog_data%var(iLookPROG%mLayerVolFracLiq)%dat, & ! volumetric fraction of liquid water in the snow+soil domain (-)
mLayerDepth => prog_data%var(iLookPROG%mLayerDepth)%dat & ! depth of each snow+soil layer (m)
- ) ! associations to variables in data structures
+ ) ! associations to variables in data structures
- ! * compute change in canopy ice content due to sublimation...
- ! ------------------------------------------------------------
- if(computeVegFlux)then
-
- ! remove mass of ice on the canopy
- scalarCanopyIce = scalarCanopyIce + scalarCanopySublimation*dt_sub
+ ! * compute change in canopy ice content due to sublimation...
+ ! ------------------------------------------------------------
+ if(computeVegFlux)then
- ! if removed all ice, take the remaining sublimation from water
- if(scalarCanopyIce < 0._dp)then
- scalarCanopyLiq = scalarCanopyLiq + scalarCanopyIce
- scalarCanopyIce = 0._dp
- endif
+ ! remove mass of ice on the canopy
+ scalarCanopyIce = scalarCanopyIce + scalarCanopySublimation*dt_sub
- ! 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
+ ! if removed all ice, take the remaining sublimation from water
+ if(scalarCanopyIce < 0._dp)then
+ scalarCanopyLiq = scalarCanopyLiq + scalarCanopyIce
+ scalarCanopyIce = 0._dp
+ endif
- end if ! (if computing the vegetation flux)
+ ! 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)
- !********** 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
+ 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; message="computSnowDepth";return; end if
! process the flag for too much sublimation
if(tooMuchSublim)then
@@ -983,161 +976,42 @@ subroutine coupled_em(&
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
+ err=20; return
endif
-
! try again
cycle substeps
endif
- ! end associate sublime
+ 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
+ ! 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
- 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)
- 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 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
- end if ! sundials
- end associate sublime
+ ! 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 MAIN SOLVER ********************************************************************************
@@ -1155,7 +1029,7 @@ subroutine coupled_em(&
! check that we have completed the sub-step
if(dt_solv >= data_step-verySmall) then
- exit substeps
+ exit substeps
endif
! adjust length of the sub-step (make sure that we don't exceed the step)
@@ -1227,202 +1101,199 @@ subroutine coupled_em(&
flux_data%var(averageFlux_meta(iVar)%ixParent)%dat(:) = flux_mean%var(iVar)%dat(:)
end do
- ! ***********************************************************************************************************************************
- ! ***********************************************************************************************************************************
- ! ***********************************************************************************************************************************
- ! ***********************************************************************************************************************************
-
- ! ---
- ! *** balance checks...
- ! ---------------------
-
- ! save the average compression and melt pond storage in the data structures
- prog_data%var(iLookPROG%scalarSfcMeltPond)%dat(1) = sfcMeltPond
-
- ! associate local variables with information in the data structures
- associate(&
- ! model forcing
- scalarSnowfall => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarSnowfall) )%dat(1) ,& ! computed snowfall rate (kg m-2 s-1)
- scalarRainfall => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarRainfall) )%dat(1) ,& ! computed rainfall rate (kg m-2 s-1)
- ! canopy fluxes
- averageThroughfallSnow => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarThroughfallSnow) )%dat(1) ,& ! snow that reaches the ground without ever touching the canopy (kg m-2 s-1)
- averageThroughfallRain => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarThroughfallRain) )%dat(1) ,& ! rain that reaches the ground without ever touching the canopy (kg m-2 s-1)
- averageCanopySnowUnloading => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarCanopySnowUnloading))%dat(1) ,& ! unloading of snow from the vegetion canopy (kg m-2 s-1)
- averageCanopyLiqDrainage => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarCanopyLiqDrainage) )%dat(1) ,& ! drainage of liquid water from the vegetation canopy (kg m-2 s-1)
- averageCanopySublimation => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarCanopySublimation) )%dat(1) ,& ! canopy sublimation/frost (kg m-2 s-1)
- averageCanopyEvaporation => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarCanopyEvaporation) )%dat(1) ,& ! canopy evaporation/condensation (kg m-2 s-1)
- ! snow fluxes
- averageSnowSublimation => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarSnowSublimation) )%dat(1) ,& ! sublimation from the snow surface (kg m-2 s-1)
- averageSnowDrainage => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarSnowDrainage) )%dat(1) ,& ! drainage from the bottom of the snowpack (m s-1)
- ! soil fluxes
- averageSoilInflux => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarInfiltration) )%dat(1) ,& ! influx of water at the top of the soil profile (m s-1)
- averageSoilDrainage => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarSoilDrainage) )%dat(1) ,& ! drainage from the bottom of the soil profile (m s-1)
- averageSoilBaseflow => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarSoilBaseflow) )%dat(1) ,& ! total baseflow from throughout the soil profile (m s-1)
- averageGroundEvaporation => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarGroundEvaporation) )%dat(1) ,& ! soil evaporation (kg m-2 s-1)
- averageCanopyTranspiration => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarCanopyTranspiration))%dat(1) ,& ! canopy transpiration (kg m-2 s-1)
- ! 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)
- ! error tolerance
- absConvTol_liquid => mpar_data%var(iLookPARAM%absConvTol_liquid)%dat(1) ,& ! absolute convergence tolerance for vol frac liq water (-)
- 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)
- ) ! (association of local variables with information in the data structures
-
- ! -----
- ! * balance checks for the canopy...
- ! ----------------------------------
-
- ! if computing the vegetation flux
- if(computeVegFlux)then
-
- ! canopy water balance
- balanceCanopyWater1 = scalarCanopyLiq + scalarCanopyIce
-
- ! balance checks for the canopy
- ! NOTE: need to put the balance checks in the sub-step loop so that we can re-compute if necessary
- scalarCanopyWatBalError = balanceCanopyWater1 - (balanceCanopyWater0 + (scalarSnowfall - averageThroughfallSnow)*data_step + (scalarRainfall - averageThroughfallRain)*data_step &
- - averageCanopySnowUnloading*data_step - averageCanopyLiqDrainage*data_step + averageCanopySublimation*data_step + averageCanopyEvaporation*data_step)
- if(abs(scalarCanopyWatBalError) > absConvTol_liquid*iden_water*10._dp)then
- print*, '** canopy water balance error:'
- write(*,'(a,1x,f20.10)') 'data_step = ', data_step
- write(*,'(a,1x,f20.10)') 'balanceCanopyWater0 = ', balanceCanopyWater0
- write(*,'(a,1x,f20.10)') 'balanceCanopyWater1 = ', balanceCanopyWater1
- write(*,'(a,1x,f20.10)') 'scalarSnowfall = ', scalarSnowfall
- write(*,'(a,1x,f20.10)') 'scalarRainfall = ', scalarRainfall
- write(*,'(a,1x,f20.10)') '(scalarSnowfall - averageThroughfallSnow) = ', (scalarSnowfall - averageThroughfallSnow)!*data_step
- write(*,'(a,1x,f20.10)') '(scalarRainfall - averageThroughfallRain) = ', (scalarRainfall - averageThroughfallRain)!*data_step
- write(*,'(a,1x,f20.10)') 'averageCanopySnowUnloading = ', averageCanopySnowUnloading!*data_step
- write(*,'(a,1x,f20.10)') 'averageCanopyLiqDrainage = ', averageCanopyLiqDrainage!*data_step
- write(*,'(a,1x,f20.10)') 'averageCanopySublimation = ', averageCanopySublimation!*data_step
- write(*,'(a,1x,f20.10)') 'averageCanopyEvaporation = ', averageCanopyEvaporation!*data_step
- write(*,'(a,1x,f20.10)') 'scalarCanopyWatBalError = ', scalarCanopyWatBalError
- message=trim(message)//'canopy hydrology does not balance'
- print*, message
- err=20; return
- end if
-
- endif ! if computing the vegetation flux
-
- ! -----
- ! * balance checks for SWE...
- ! ---------------------------
-
- ! recompute snow depth (m) and SWE (kg m-2)
- 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
-
- ! check the individual layers
- if(printBalance .and. nSnow>0)then
- write(*,'(a,1x,10(f12.8,1x))') 'liqSnowInit = ', liqSnowInit
- write(*,'(a,1x,10(f12.8,1x))') 'volFracLiq = ', prog_data%var(iLookPROG%mLayerVolFracLiq)%dat(1:nSnow)
- write(*,'(a,1x,10(f12.8,1x))') 'iLayerLiqFluxSnow = ', flux_data%var(iLookFLUX%iLayerLiqFluxSnow)%dat*iden_water*data_step
- write(*,'(a,1x,10(f12.8,1x))') 'mLayerLiqFluxSnow = ', flux_data%var(iLookFLUX%mLayerLiqFluxSnow)%dat*data_step
- write(*,'(a,1x,10(f12.8,1x))') 'change volFracLiq = ', prog_data%var(iLookPROG%mLayerVolFracLiq)%dat(1:nSnow) - liqSnowInit
- deallocate(liqSnowInit, stat=err)
- if(err/=0)then
- message=trim(message)//'unable to deallocate space for the initial volumetric liquid water content of snow'
- print*, message
- err=20; return
- endif
- endif
-
- ! check SWE
- if(nSnow>0)then
- effSnowfall = averageThroughfallSnow + averageCanopySnowUnloading
- effRainfall = averageThroughfallRain + averageCanopyLiqDrainage
- newSWE = prog_data%var(iLookPROG%scalarSWE)%dat(1)
- delSWE = newSWE - (oldSWE - sfcMeltPond)
- massBalance = delSWE - (effSnowfall + effRainfall + averageSnowSublimation - averageSnowDrainage*iden_water)*data_step
-
- if(abs(massBalance) > absConvTol_liquid*iden_water*10._dp)then
- print*, 'nSnow = ', nSnow
- print*, 'nSub = ', nSub
- write(*,'(a,1x,f20.10)') 'data_step = ', data_step
- write(*,'(a,1x,f20.10)') 'oldSWE = ', oldSWE
- write(*,'(a,1x,f20.10)') 'newSWE = ', newSWE
- write(*,'(a,1x,f20.10)') 'delSWE = ', delSWE
- write(*,'(a,1x,f20.10)') 'effRainfall = ', effRainfall*data_step
- write(*,'(a,1x,f20.10)') 'effSnowfall = ', effSnowfall*data_step
- write(*,'(a,1x,f20.10)') 'sublimation = ', averageSnowSublimation*data_step
- write(*,'(a,1x,f20.10)') 'snwDrainage = ', averageSnowDrainage*iden_water*data_step
- write(*,'(a,1x,f20.10)') 'sfcMeltPond = ', sfcMeltPond
- write(*,'(a,1x,f20.10)') 'massBalance = ', massBalance
- message=trim(message)//'SWE does not balance'
- print*,message
- err=20; return
- endif ! if failed mass balance check
- endif ! if snow layers exist
-
- ! -----
- ! * balance checks for soil...
- ! ----------------------------
+ ! ***********************************************************************************************************************************
+ ! ***********************************************************************************************************************************
+ ! ***********************************************************************************************************************************
+ ! ***********************************************************************************************************************************
+
+ ! ---
+ ! *** balance checks...
+ ! ---------------------
+
+ ! save the average compression and melt pond storage in the data structures
+ prog_data%var(iLookPROG%scalarSfcMeltPond)%dat(1) = sfcMeltPond
+
+ ! associate local variables with information in the data structures
+ associate(&
+ ! model forcing
+ scalarSnowfall => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarSnowfall) )%dat(1) ,& ! computed snowfall rate (kg m-2 s-1)
+ scalarRainfall => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarRainfall) )%dat(1) ,& ! computed rainfall rate (kg m-2 s-1)
+ ! canopy fluxes
+ averageThroughfallSnow => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarThroughfallSnow) )%dat(1) ,& ! snow that reaches the ground without ever touching the canopy (kg m-2 s-1)
+ averageThroughfallRain => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarThroughfallRain) )%dat(1) ,& ! rain that reaches the ground without ever touching the canopy (kg m-2 s-1)
+ averageCanopySnowUnloading => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarCanopySnowUnloading))%dat(1) ,& ! unloading of snow from the vegetion canopy (kg m-2 s-1)
+ averageCanopyLiqDrainage => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarCanopyLiqDrainage) )%dat(1) ,& ! drainage of liquid water from the vegetation canopy (kg m-2 s-1)
+ averageCanopySublimation => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarCanopySublimation) )%dat(1) ,& ! canopy sublimation/frost (kg m-2 s-1)
+ averageCanopyEvaporation => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarCanopyEvaporation) )%dat(1) ,& ! canopy evaporation/condensation (kg m-2 s-1)
+ ! snow fluxes
+ averageSnowSublimation => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarSnowSublimation) )%dat(1) ,& ! sublimation from the snow surface (kg m-2 s-1)
+ averageSnowDrainage => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarSnowDrainage) )%dat(1) ,& ! drainage from the bottom of the snowpack (m s-1)
+ ! soil fluxes
+ averageSoilInflux => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarInfiltration) )%dat(1) ,& ! influx of water at the top of the soil profile (m s-1)
+ averageSoilDrainage => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarSoilDrainage) )%dat(1) ,& ! drainage from the bottom of the soil profile (m s-1)
+ averageSoilBaseflow => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarSoilBaseflow) )%dat(1) ,& ! total baseflow from throughout the soil profile (m s-1)
+ averageGroundEvaporation => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarGroundEvaporation) )%dat(1) ,& ! soil evaporation (kg m-2 s-1)
+ averageCanopyTranspiration => flux_mean%var(childFLUX_MEAN(iLookFLUX%scalarCanopyTranspiration))%dat(1) ,& ! canopy transpiration (kg m-2 s-1)
+ ! 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)
+ ! error tolerance
+ absConvTol_liquid => mpar_data%var(iLookPARAM%absConvTol_liquid)%dat(1) ,& ! absolute convergence tolerance for vol frac liq water (-)
+ 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)
+ ) ! (association of local variables with information in the data structures
+
+ ! -----
+ ! * balance checks for the canopy...
+ ! ----------------------------------
+
+ if (model_decisions(iLookDECISIONS%num_method)%iDecision==bEuler) checkMassBalance = .true. ! convergence criteria for bEuler
+ if (model_decisions(iLookDECISIONS%num_method)%iDecision==sundials) checkMassBalance = .false. ! sundials does not use this criteria
+
+ ! if computing the vegetation flux
+ if(computeVegFlux)then
+
+ ! canopy water balance
+ balanceCanopyWater1 = scalarCanopyLiq + scalarCanopyIce
+
+ ! balance checks for the canopy
+ ! NOTE: need to put the balance checks in the sub-step loop so that we can re-compute if necessary
+ scalarCanopyWatBalError = balanceCanopyWater1 - (balanceCanopyWater0 + (scalarSnowfall - averageThroughfallSnow)*data_step + (scalarRainfall - averageThroughfallRain)*data_step &
+ - averageCanopySnowUnloading*data_step - averageCanopyLiqDrainage*data_step + averageCanopySublimation*data_step + averageCanopyEvaporation*data_step)
+ if(abs(scalarCanopyWatBalError) > absConvTol_liquid*iden_water*10._rkind .and. checkMassBalance)then
+ print*, '** canopy water balance error:'
+ write(*,'(a,1x,f20.10)') 'data_step = ', data_step
+ write(*,'(a,1x,f20.10)') 'balanceCanopyWater0 = ', balanceCanopyWater0
+ write(*,'(a,1x,f20.10)') 'balanceCanopyWater1 = ', balanceCanopyWater1
+ write(*,'(a,1x,f20.10)') 'scalarSnowfall = ', scalarSnowfall
+ write(*,'(a,1x,f20.10)') 'scalarRainfall = ', scalarRainfall
+ write(*,'(a,1x,f20.10)') '(scalarSnowfall - averageThroughfallSnow) = ', (scalarSnowfall - averageThroughfallSnow)!*data_step
+ write(*,'(a,1x,f20.10)') '(scalarRainfall - averageThroughfallRain) = ', (scalarRainfall - averageThroughfallRain)!*data_step
+ write(*,'(a,1x,f20.10)') 'averageCanopySnowUnloading = ', averageCanopySnowUnloading!*data_step
+ write(*,'(a,1x,f20.10)') 'averageCanopyLiqDrainage = ', averageCanopyLiqDrainage!*data_step
+ write(*,'(a,1x,f20.10)') 'averageCanopySublimation = ', averageCanopySublimation!*data_step
+ write(*,'(a,1x,f20.10)') 'averageCanopyEvaporation = ', averageCanopyEvaporation!*data_step
+ write(*,'(a,1x,f20.10)') 'scalarCanopyWatBalError = ', scalarCanopyWatBalError
+ message=trim(message)//'canopy hydrology does not balance'
+ err=20; return
+ end if
- ! compute the liquid water and ice content at the end of the time step
- scalarTotalSoilLiq = sum(iden_water*mLayerVolFracLiq(1:nSoil)*mLayerDepth(1:nSoil))
- scalarTotalSoilIce = sum(iden_water*mLayerVolFracIce(1:nSoil)*mLayerDepth(1:nSoil)) ! NOTE: no expansion of soil, hence use iden_water
+ endif ! if computing the vegetation flux
- ! get the total water in the soil (liquid plus ice) at the end of the time step (kg m-2)
- balanceSoilWater1 = scalarTotalSoilLiq + scalarTotalSoilIce
+ ! -----
+ ! * balance checks for SWE...
+ ! ---------------------------
- ! get the total aquifer storage at the start of the time step (kg m-2)
- balanceAquifer1 = scalarAquiferStorage*iden_water
+ ! recompute snow depth (m) and SWE (kg m-2)
+ 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
- ! get the input and output to/from the soil zone (kg m-2)
- balanceSoilInflux = averageSoilInflux*iden_water*data_step
- balanceSoilBaseflow = averageSoilBaseflow*iden_water*data_step
- balanceSoilDrainage = averageSoilDrainage*iden_water*data_step
- balanceSoilET = (averageCanopyTranspiration + averageGroundEvaporation)*data_step
+ ! check the individual layers
+ if(printBalance .and. nSnow>0)then
+ write(*,'(a,1x,10(f12.8,1x))') 'liqSnowInit = ', liqSnowInit
+ write(*,'(a,1x,10(f12.8,1x))') 'volFracLiq = ', prog_data%var(iLookPROG%mLayerVolFracLiq)%dat(1:nSnow)
+ write(*,'(a,1x,10(f12.8,1x))') 'iLayerLiqFluxSnow = ', flux_data%var(iLookFLUX%iLayerLiqFluxSnow)%dat*iden_water*data_step
+ write(*,'(a,1x,10(f12.8,1x))') 'mLayerLiqFluxSnow = ', flux_data%var(iLookFLUX%mLayerLiqFluxSnow)%dat*data_step
+ write(*,'(a,1x,10(f12.8,1x))') 'change volFracLiq = ', prog_data%var(iLookPROG%mLayerVolFracLiq)%dat(1:nSnow) - liqSnowInit
+ deallocate(liqSnowInit, stat=err)
+ if(err/=0)then
+ message=trim(message)//'unable to deallocate space for the initial volumetric liquid water content of snow'
+ err=20; return
+ endif
+ endif
- ! check the individual layers
- if(printBalance)then
- write(*,'(a,1x,10(f12.8,1x))') 'liqSoilInit = ', liqSoilInit
- write(*,'(a,1x,10(f12.8,1x))') 'volFracLiq = ', mLayerVolFracLiq
- write(*,'(a,1x,10(f12.8,1x))') 'iLayerLiqFluxSoil = ', flux_data%var(iLookFLUX%iLayerLiqFluxSoil)%dat*iden_water*data_step
- write(*,'(a,1x,10(f12.8,1x))') 'mLayerLiqFluxSoil = ', flux_data%var(iLookFLUX%mLayerLiqFluxSoil)%dat*data_step
- write(*,'(a,1x,10(f12.8,1x))') 'change volFracLiq = ', mLayerVolFracLiq - liqSoilInit
- deallocate(liqSoilInit, stat=err)
- if(err/=0)then
- message=trim(message)//'unable to deallocate space for the initial soil moisture'
- err=20; return
- print*, message
- endif
- endif
+ ! check SWE
+ if(nSnow>0)then
+ effSnowfall = averageThroughfallSnow + averageCanopySnowUnloading
+ effRainfall = averageThroughfallRain + averageCanopyLiqDrainage
+ newSWE = prog_data%var(iLookPROG%scalarSWE)%dat(1)
+ delSWE = newSWE - (oldSWE - sfcMeltPond)
+ massBalance = delSWE - (effSnowfall + effRainfall + averageSnowSublimation - averageSnowDrainage*iden_water)*data_step
+ if(abs(massBalance) > absConvTol_liquid*iden_water*10._rkind .and. checkMassBalance)then
+ print*, 'nSnow = ', nSnow
+ print*, 'nSub = ', nSub
+ write(*,'(a,1x,f20.10)') 'data_step = ', data_step
+ write(*,'(a,1x,f20.10)') 'oldSWE = ', oldSWE
+ write(*,'(a,1x,f20.10)') 'newSWE = ', newSWE
+ write(*,'(a,1x,f20.10)') 'delSWE = ', delSWE
+ write(*,'(a,1x,f20.10)') 'effRainfall = ', effRainfall*data_step
+ write(*,'(a,1x,f20.10)') 'effSnowfall = ', effSnowfall*data_step
+ write(*,'(a,1x,f20.10)') 'sublimation = ', averageSnowSublimation*data_step
+ write(*,'(a,1x,f20.10)') 'snwDrainage = ', averageSnowDrainage*iden_water*data_step
+ write(*,'(a,1x,f20.10)') 'sfcMeltPond = ', sfcMeltPond
+ write(*,'(a,1x,f20.10)') 'massBalance = ', massBalance
+ message=trim(message)//'SWE does not balance'
+ err=20; return
+ endif ! if failed mass balance check
+ endif ! if snow layers exist
+
+ ! -----
+ ! * balance checks for soil...
+ ! ----------------------------
+
+ ! compute the liquid water and ice content at the end of the time step
+ scalarTotalSoilLiq = sum(iden_water*mLayerVolFracLiq(1:nSoil)*mLayerDepth(1:nSoil))
+ scalarTotalSoilIce = sum(iden_water*mLayerVolFracIce(1:nSoil)*mLayerDepth(1:nSoil)) ! NOTE: no expansion of soil, hence use iden_water
+
+ ! get the total water in the soil (liquid plus ice) at the end of the time step (kg m-2)
+ balanceSoilWater1 = scalarTotalSoilLiq + scalarTotalSoilIce
+
+ ! get the total aquifer storage at the start of the time step (kg m-2)
+ balanceAquifer1 = scalarAquiferStorage*iden_water
+
+ ! get the input and output to/from the soil zone (kg m-2)
+ balanceSoilInflux = averageSoilInflux*iden_water*data_step
+ balanceSoilBaseflow = averageSoilBaseflow*iden_water*data_step
+ balanceSoilDrainage = averageSoilDrainage*iden_water*data_step
+ balanceSoilET = (averageCanopyTranspiration + averageGroundEvaporation)*data_step
+
+ ! check the individual layers
+ if(printBalance)then
+ write(*,'(a,1x,10(f12.8,1x))') 'liqSoilInit = ', liqSoilInit
+ write(*,'(a,1x,10(f12.8,1x))') 'volFracLiq = ', mLayerVolFracLiq
+ write(*,'(a,1x,10(f12.8,1x))') 'iLayerLiqFluxSoil = ', flux_data%var(iLookFLUX%iLayerLiqFluxSoil)%dat*iden_water*data_step
+ write(*,'(a,1x,10(f12.8,1x))') 'mLayerLiqFluxSoil = ', flux_data%var(iLookFLUX%mLayerLiqFluxSoil)%dat*data_step
+ write(*,'(a,1x,10(f12.8,1x))') 'change volFracLiq = ', mLayerVolFracLiq - liqSoilInit
+ deallocate(liqSoilInit, stat=err)
+ if(err/=0)then
+ message=trim(message)//'unable to deallocate space for the initial soil moisture'
+ err=20; return
+ endif
+ endif
- ! check the soil water balance
- scalarSoilWatBalError = balanceSoilWater1 - (balanceSoilWater0 + (balanceSoilInflux + balanceSoilET - balanceSoilBaseflow - balanceSoilDrainage - totalSoilCompress) )
- if(abs(scalarSoilWatBalError) > absConvTol_liquid*iden_water*10._dp)then ! NOTE: kg m-2, so need coarse tolerance to account for precision issues
- write(*,*) 'solution method = ', ixSolution
- write(*,'(a,1x,f20.10)') 'data_step = ', data_step
- write(*,'(a,1x,f20.10)') 'totalSoilCompress = ', totalSoilCompress
- write(*,'(a,1x,f20.10)') 'scalarTotalSoilLiq = ', scalarTotalSoilLiq
- write(*,'(a,1x,f20.10)') 'scalarTotalSoilIce = ', scalarTotalSoilIce
- write(*,'(a,1x,f20.10)') 'balanceSoilWater0 = ', balanceSoilWater0
- write(*,'(a,1x,f20.10)') 'balanceSoilWater1 = ', balanceSoilWater1
- write(*,'(a,1x,f20.10)') 'balanceSoilInflux = ', balanceSoilInflux
- write(*,'(a,1x,f20.10)') 'balanceSoilBaseflow = ', balanceSoilBaseflow
- write(*,'(a,1x,f20.10)') 'balanceSoilDrainage = ', balanceSoilDrainage
- write(*,'(a,1x,f20.10)') 'balanceSoilET = ', balanceSoilET
- write(*,'(a,1x,f20.10)') 'scalarSoilWatBalError = ', scalarSoilWatBalError
- write(*,'(a,1x,f20.10)') 'scalarSoilWatBalError = ', scalarSoilWatBalError/iden_water
- write(*,'(a,1x,f20.10)') 'absConvTol_liquid = ', absConvTol_liquid
- ! error control
- message=trim(message)//'soil hydrology does not balance'
- print*, message
- err=20; return
- end if
+ ! check the soil water balance
+ scalarSoilWatBalError = balanceSoilWater1 - (balanceSoilWater0 + (balanceSoilInflux + balanceSoilET - balanceSoilBaseflow - balanceSoilDrainage - totalSoilCompress) )
+ if(abs(scalarSoilWatBalError) > absConvTol_liquid*iden_water*10._rkind .and. checkMassBalance)then ! NOTE: kg m-2, so need coarse tolerance to account for precision issues
+ write(*,*) 'solution method = ', ixSolution
+ write(*,'(a,1x,f20.10)') 'data_step = ', data_step
+ write(*,'(a,1x,f20.10)') 'totalSoilCompress = ', totalSoilCompress
+ write(*,'(a,1x,f20.10)') 'scalarTotalSoilLiq = ', scalarTotalSoilLiq
+ write(*,'(a,1x,f20.10)') 'scalarTotalSoilIce = ', scalarTotalSoilIce
+ write(*,'(a,1x,f20.10)') 'balanceSoilWater0 = ', balanceSoilWater0
+ write(*,'(a,1x,f20.10)') 'balanceSoilWater1 = ', balanceSoilWater1
+ write(*,'(a,1x,f20.10)') 'balanceSoilInflux = ', balanceSoilInflux
+ write(*,'(a,1x,f20.10)') 'balanceSoilBaseflow = ', balanceSoilBaseflow
+ write(*,'(a,1x,f20.10)') 'balanceSoilDrainage = ', balanceSoilDrainage
+ write(*,'(a,1x,f20.10)') 'balanceSoilET = ', balanceSoilET
+ write(*,'(a,1x,f20.10)') 'scalarSoilWatBalError = ', scalarSoilWatBalError
+ write(*,'(a,1x,f20.10)') 'scalarSoilWatBalError = ', scalarSoilWatBalError/iden_water
+ write(*,'(a,1x,f20.10)') 'absConvTol_liquid = ', absConvTol_liquid
+ ! error control
+ message=trim(message)//'soil hydrology does not balance'
+ err=20; return
+ end if
- ! end association of local variables with information in the data structures
- end associate
+ ! end association of local variables with information in the data structures
+ end associate
! end association to canopy depth
end associate canopy
@@ -1440,13 +1311,12 @@ subroutine coupled_em(&
end if
iLayer = nSnow+1
- !print*, 'nsub, mLayerTemp(iLayer), mLayerVolFracIce(iLayer) = ', nsub, mLayerTemp(iLayer), mLayerVolFracIce(iLayer)
- !print*, 'nsub = ', nsub
if(nsub>50000)then
write(message,'(a,i0)') trim(cmessage)//'number of sub-steps > 50000 for HRU ', indxHRU
err=20; return
end if
+
end subroutine coupled_em