From 8f945222ec3ee5fbd5cd24ff1ea5fd6745e78cba Mon Sep 17 00:00:00 2001
From: Kyle <kyle.c.klenk@gmail.com>
Date: Mon, 21 Nov 2022 21:25:52 +0000
Subject: [PATCH] changed summaActors_fileManager to summaFileManager
This allows the easier hookup of the old summa which has the ladder module name
---
build/makefile_old_summa | 435 +++++
build/makefile_sundials | 2 +-
.../fortran_code/cppwrap_fileAccess.f90 | 6 +-
.../fortran_code/read_attribute.f90 | 8 +-
.../fortran_code/read_forcing.f90 | 4 +-
.../fortran_code/read_param.f90 | 4 +-
.../hru_actor/fortran_code/hru_actor.f90 | 2 +-
build/source/driver/summaActors_restart.f90 | 6 +-
build/source/dshare/popMetadat.f90 | 4 +-
build/source/engine/coupled_em.f90 | 223 ++-
build/source/engine/derivforce.f90 | 2 +-
build/source/engine/ffile_info.f90 | 8 +-
build/source/engine/mDecisions.f90 | 6 +-
build/source/engine/read_dimension.f90 | 4 +-
build/source/engine/read_pinit.f90 | 2 +-
build/source/engine/sundials/coupled_em.f90 | 1516 +++++++++++++++++
build/source/engine/vegPhenlgy.f90 | 2 -
.../source/hookup/summaActors_FileManager.f90 | 4 +-
build/source/netcdf/def_output.f90 | 2 +-
build/source/netcdf/read_icondActors.f90 | 6 +-
build/summa | 1 +
21 files changed, 2098 insertions(+), 149 deletions(-)
create mode 100644 build/makefile_old_summa
mode change 100755 => 100644 build/source/engine/coupled_em.f90
create mode 100755 build/source/engine/sundials/coupled_em.f90
create mode 160000 build/summa
diff --git a/build/makefile_old_summa b/build/makefile_old_summa
new file mode 100644
index 0000000..b9f5c50
--- /dev/null
+++ b/build/makefile_old_summa
@@ -0,0 +1,435 @@
+#### parent directory of the 'build' directory ####
+ROOT_DIR = /Summa-Actors
+
+#### Compilers ####
+FC = gfortran # Fortran
+CC = g++ # C++
+
+DIR_SUNDIALS=/code/sundials/instdir
+INC_SUNDIALS=-I$(DIR_SUNDIALS)/include -I$(DIR_SUNDIALS)/fortran
+LIB_SUNDIALS=-L$(DIR_SUNDIALS)/lib -lsundials_fnvecmanyvector_mod -lsundials_fida_mod -lsundials_fnvecserial_mod -lsundials_fsunlinsoldense_mod -lsundials_fsunmatrixdense_mod
+
+
+#### Includes AND Libraries ####
+INCLUDES = -I/usr/include -I/usr/local/include $(INC_SUNDIALS)
+LIBRARIES = -L/usr/lib -L/usr/local/lib -lnetcdff -lopenblas $(LIB_SUNDIALS)
+
+ACTORS_INCLUDES = -I/usr/include -I/usr/local/include $(INC_SUNDIALS)
+ACTORS_LIBRARIES = -L/usr/lib -L/usr/local/lib -L/Summa-Actors/bin -lcaf_core -lcaf_io -lsumma -lopenblas -lnetcdff $(LIB_SUNDIALS)
+
+
+# Production runs
+FLAGS_NOAH = -g -O0 -ffree-form -ffree-line-length-none -fmax-errors=0 -fPIC -Wfatal-errors
+FLAGS_COMM = -g -O0 -ffree-line-length-none -fmax-errors=0 -fPIC -Wfatal-errors
+FLAGS_SUMMA = -g -O0 -ffree-line-length-none -fmax-errors=0 -fPIC -Wfatal-errors
+FLAGS_ACTORS = -g -O0 -Wfatal-errors -std=c++17
+
+# Debug runs
+# FLAGS_NOAH = -g -O0 -ffree-form -ffree-line-length-none -fmax-errors=0 -fbacktrace -Wno-unused -Wno-unused-dummy-argument -fPIC
+# FLAGS_COMM = -g -O0 -Wall -ffree-line-length-none -fmax-errors=0 -fbacktrace -fcheck=bounds -fPIC
+# FLAGS_SUMMA = -g -O0 -Wall -ffree-line-length-none -fmax-errors=0 -fbacktrace -fcheck=bounds -fPIC
+# FLAGS_ACTORS = -g -O0 -Wall -std=c++17
+
+
+
+#========================================================================
+# PART 1: Define directory paths
+#========================================================================
+
+# Core directory that contains source code
+F_KORE_DIR = $(ROOT_DIR)/build/source
+
+# Location of the compiled modules
+MOD_PATH = $(ROOT_DIR)/build
+
+# Define the directory for the executables
+EXE_PATH = $(ROOT_DIR)/bin
+
+####################################################################################################
+###################################### Assemble Fortran Files ######################################
+####################################################################################################
+# Define directories
+DRIVER_DIR = $(F_KORE_DIR)/driver
+HOOKUP_DIR = $(F_KORE_DIR)/hookup
+NETCDF_DIR = $(F_KORE_DIR)/netcdf
+DSHARE_DIR = $(F_KORE_DIR)/dshare
+NUMREC_DIR = $(F_KORE_DIR)/numrec
+NOAHMP_DIR = $(F_KORE_DIR)/noah-mp
+ENGINE_DIR = $(F_KORE_DIR)/engine
+ACTORS_DIR = $(F_KORE_DIR)/actors
+
+SUMMA_ENGINE_DIR = $(MOD_PATH)/summa/build/source/engine
+
+
+JOB_ACTOR_DIR = $(ACTORS_DIR)/job_actor
+FILE_ACCESS_DIR = $(ACTORS_DIR)/file_access_actor/fortran_code
+HRU_ACTOR_DIR = $(ACTORS_DIR)/hru_actor/fortran_code
+GRU_ACTOR_DIR = $(ACTORS_DIR)/gru_actor
+
+# utilities
+SUMMA_NRUTIL= \
+ nrtype.f90 \
+ f2008funcs.f90 \
+ nr_utility.f90 \
+
+NRUTIL = $(patsubst %, $(ENGINE_DIR)/%, $(SUMMA_NRUTIL))
+
+# Numerical recipes procedures
+# NOTE: all numerical recipes procedures are now replaced with free versions
+SUMMA_NRPROC= \
+ expIntegral.f90 \
+ spline_int.f90
+NRPROC = $(patsubst %, $(ENGINE_DIR)/%, $(SUMMA_NRPROC))
+
+# Hook-up modules (set files and directory paths)
+SUMMA_HOOKUP= \
+ ascii_util.f90 \
+ summaActors_FileManager.f90
+HOOKUP = $(patsubst %, $(HOOKUP_DIR)/%, $(SUMMA_HOOKUP))
+
+# Data modules
+SUMMA_DATAMS= \
+ multiconst.f90 \
+ var_lookup.f90 \
+ data_types.f90 \
+ globalData.f90 \
+ flxMapping.f90 \
+ get_ixname.f90
+DATAMS = $(patsubst %, $(DSHARE_DIR)/%, $(SUMMA_DATAMS))
+
+SUMMA_DEPEND_ON_FILEMANAGER= \
+ popMetadat.f90 \
+ outpt_stat.f90
+DEPEND_ON_FILEMANAGER = $(patsubst %, $(DSHARE_DIR)/%, $(SUMMA_DEPEND_ON_FILEMANAGER))
+
+
+
+
+# utility modules
+TIME_UTILS = \
+ time_utils.f90
+
+TIME_UTILS_UTILMS = $(patsubst %, $(ENGINE_DIR)/%, $(TIME_UTILS))
+
+M_DECISIONS = \
+ mDecisions.f90 \
+
+M_DECISIONS_UTILMS = $(patsubst %, $(SUMMA_ENGINE_DIR)/%, $(M_DECISIONS))
+
+
+SUMMA_UTILMS= \
+ snow_utils.f90 \
+ soil_utils.f90 \
+ sundials/soil_utilsAddSundials.f90 \
+ updatState.f90 \
+ sundials/updatStateSundials.f90 \
+ matrixOper.f90
+UTILMS = $(patsubst %, $(ENGINE_DIR)/%, $(SUMMA_UTILMS))
+
+# Model guts
+SUMMA_MODGUT= \
+ MODGUT = $(patsubst %, $(ENGINE_DIR)/%, $(SUMMA_MODGUT))
+
+# Solver - Split up to compiling with old summa files vs new summa files
+
+VEG_PHENLGY= \
+ vegPhenlgy.f90
+
+VEG_PHENLGY_SOLVER = $(patsubst %, $(ENGINE_DIR)/%, $(VEG_PHENLGY))
+
+
+SUMMA_SOLVER= \
+ diagn_evar.f90 \
+ stomResist.f90 \
+ groundwatr.f90 \
+ vegSWavRad.f90 \
+ vegNrgFlux.f90 \
+ ssdNrgFlux.f90 \
+ vegLiqFlux.f90 \
+ snowLiqFlx.f90 \
+ soilLiqFlx.f90 \
+ bigAquifer.f90 \
+ computFlux.f90 \
+ computResid.f90 \
+ computJacob.f90 \
+ eval8summa.f90 \
+ summaSolve.f90 \
+ systemSolv.f90 \
+ opSplittin.f90
+
+SOLVER = $(patsubst %, $(SUMMA_ENGINE_DIR)/%, $(SUMMA_SOLVER))
+
+SUMMA_ACTORS_SOLVER = \
+ coupled_em.f90
+
+COUPLED_EM = $(patsubst %, $(ENGINE_DIR)/%, $(SUMMA_ACTORS_SOLVER))
+
+
+
+# Interface code for Fortran-C++
+SUMMA_INTERFACE= \
+ cppwrap_datatypes.f90 \
+ cppwrap_auxiliary.f90 \
+ cppwrap_metadata.f90 \
+
+INTERFACE = $(patsubst %, $(ACTORS_DIR)/global/%, $(SUMMA_INTERFACE))
+
+
+SUMMA_FILEACCESS_INTERFACE = \
+ cppwrap_fileAccess.f90 \
+ read_attribute.f90 \
+ read_forcing.f90 \
+ read_param.f90 \
+ write_to_netcdf.f90
+
+FILEACCESS_INTERFACE = $(patsubst %, $(FILE_ACCESS_DIR)/%, $(SUMMA_FILEACCESS_INTERFACE))
+
+SUMMA_JOB_INTERFACE = \
+ job_actor.f90
+
+JOB_INTERFACE = $(patsubst %, $(JOB_ACTOR_DIR)/%, $(SUMMA_JOB_INTERFACE))
+
+SUMMA_HRU_INTERFACE = \
+ cppwrap_hru.f90 \
+ hru_actor.f90 \
+ init_hru_actor.f90 \
+
+
+HRU_INTERFACE = $(patsubst %, $(HRU_ACTOR_DIR)/%, $(SUMMA_HRU_INTERFACE))
+
+SUMMA_GRU_INTERFACE = \
+ gru_actor.f90 \
+
+GRU_INTERFACE = $(patsubst %, $(GRU_ACTOR_DIR)/%, $(SUMMA_GRU_INTERFACE))
+
+
+
+# Define routines for SUMMA preliminaries
+SUMMA_PRELIM= \
+ conv_funcs.f90 \
+ sunGeomtry.f90 \
+ convE2Temp.f90 \
+ allocspaceActors.f90 \
+ checkStruc.f90 \
+ childStruc.f90 \
+ ffile_info.f90 \
+ read_dimension.f90 \
+ read_pinit.f90 \
+ pOverwrite.f90 \
+ paramCheck.f90 \
+ check_icondActors.f90 \
+ # allocspace.f90
+PRELIM = $(patsubst %, $(ENGINE_DIR)/%, $(SUMMA_PRELIM))
+
+SUMMA_NOAHMP= \
+ module_model_constants.F \
+ module_sf_noahutl.F \
+ module_sf_noahlsm.F \
+ module_sf_noahmplsm.F
+
+NOAHMP = $(patsubst %, $(NOAHMP_DIR)/%, $(SUMMA_NOAHMP))
+
+# Define routines for the SUMMA model runs
+# sundials - t2enthalpy.f90
+SUMMA_MODRUN = \
+ indexState.f90 \
+ getVectorz.f90 \
+ sundials/t2enthalpy.f90 \
+ updateVars.f90 \
+ sundials/updateVarsSundials.f90 \
+ var_derive.f90 \
+ derivforce.f90 \
+ snowAlbedo.f90 \
+ canopySnow.f90 \
+ tempAdjust.f90 \
+ snwCompact.f90 \
+ layerMerge.f90 \
+ layerDivide.f90 \
+ volicePack.f90 \
+ qTimeDelay.f90
+MODRUN = $(patsubst %, $(ENGINE_DIR)/%, $(SUMMA_MODRUN))
+
+# Define NetCDF routines
+# OutputStrucWrite is not a netcdf subroutine and should be
+# moved
+SUMMA_NETCDF = \
+ netcdf_util.f90 \
+ def_output.f90 \
+ writeOutput.f90 \
+ read_icondActors.f90
+NETCDF = $(patsubst %, $(NETCDF_DIR)/%, $(SUMMA_NETCDF))
+
+# ... stitch together common programs
+COMM_ALL = $(NRUTIL) $(NRPROC) $(DATAMS) $(INTERFACE) $(HOOKUP) $(DEPEND_ON_FILEMANAGER) \
+ $(TIME_UTILS_UTILMS) $(M_DECISIONS_UTILMS) $(UTILMS)
+# ... stitch together SUMMA programs
+SUMMA_ALL = $(NETCDF) $(PRELIM) $(MODRUN) $(VEG_PHENLGY_SOLVER) $(SOLVER) $(COUPLED_EM)
+
+# Define the driver routine
+SUMMA_DRIVER= \
+ summaActors_type.f90 \
+ summaActors_util.f90 \
+ summaActors_globalData.f90 \
+ SummaActors_setup.f90 \
+ summaActors_restart.f90 \
+ SummaActors_modelRun.f90 \
+ summaActors_alarms.f90
+
+
+DRIVER = $(patsubst %, $(DRIVER_DIR)/%, $(SUMMA_DRIVER))
+
+####################################################################################################
+###################################### Assemble Fortran Files ######################################
+####################################################################################################
+
+####################################################################################################
+######################################## Assemble C++ Files ########################################
+####################################################################################################
+
+INCLUDE_DIR = /Summa-Actors/build/includes
+SOURCE_DIR = /Summa-Actors/build/source/actors
+
+
+GLOBAL_INCLUDES = -I$(INCLUDE_DIR)/global
+GLOBAL = $(SOURCE_DIR)/global/global.cpp
+TIMEINFO = $(SOURCE_DIR)/global/timing_info.cpp
+SETTINGS_FILES = $(SOURCE_DIR)/global/settings_functions.cpp
+AUXILARY = $(SOURCE_DIR)/global/auxiliary.cpp
+
+SUMMA_ACTOR_INCLUDES = -I$(INCLUDE_DIR)/summa_actor
+SUMMA_ACTOR = $(SOURCE_DIR)/summa_actor/summa_actor.cpp
+SUMMA_CLIENT = $(SOURCE_DIR)/summa_actor/summa_client.cpp
+SUMMA_SERVER = $(SOURCE_DIR)/summa_actor/summa_server.cpp
+SUMMA_BACKUP_SERVER = $(SOURCE_DIR)/summa_actor/summa_backup_server.cpp
+
+BATCH = $(SOURCE_DIR)/summa_actor/batch/batch.cpp
+BATCH_CONTAINER = $(SOURCE_DIR)/summa_actor/batch/batch_container.cpp
+
+CLIENT = $(SOURCE_DIR)/summa_actor/client/client.cpp
+CLIENT_CONTAINER = $(SOURCE_DIR)/summa_actor/client/client_container.cpp
+
+CLIENT_BATCH = $(SOURCE_DIR)/summa_actor/batch_client.cpp
+CLIENT_BATCH_CONTAINERS = $(SOURCE_DIR)/summa_actor/batch_client_containers.cpp
+
+GRU_ACTOR_INCLUDES = -I$(INCLUDE_DIR)/gru_actor
+GRU_ACTOR = $(SOURCE_DIR)/gru_actor/gru_actor.cpp
+
+JOB_ACTOR_INCLUDES = -I$(INCLUDE_DIR)/job_actor
+JOB_ACTOR = $(SOURCE_DIR)/job_actor/job_actor.cpp
+GRUinfo = $(SOURCE_DIR)/job_actor/GRUinfo.cpp
+
+FILE_ACCESS_ACTOR_INCLUDES = -I$(INCLUDE_DIR)/file_access_actor
+FILE_ACCESS_ACTOR = $(SOURCE_DIR)/file_access_actor/cpp_code/file_access_actor.cpp
+FORCING_FILE_INFO = $(SOURCE_DIR)/file_access_actor/cpp_code/forcing_file_info.cpp
+OUTPUT_MANAGER = $(SOURCE_DIR)/file_access_actor/cpp_code/output_manager.cpp
+
+HRU_ACTOR_INCLUDES = -I$(INCLUDE_DIR)/hru_actor
+HRU_ACTOR = $(SOURCE_DIR)/hru_actor/cpp_code/hru_actor.cpp
+
+MAIN = $(F_KORE_DIR)/actors/main.cpp
+
+ACTOR_TEST = $(F_KORE_DIR)/testing/testing_main.cc
+####################################################################################################
+######################################## Assemble C++ Files ########################################
+####################################################################################################
+
+
+#========================================================================
+# PART 3: compilation
+#======================================================================
+all: fortran cpp
+
+fortran: compile_noah compile_comm compile_summa link clean_fortran
+
+cpp: compile_globals compile_hru_actor compile_gru_actor compile_file_access_actor compile_job_actor compile_summa_actor \
+ compile_summa_client compile_summa_server compile_main link_cpp clean_cpp
+
+test: actors_test actors_testLink actorsClean
+
+###################################################################################################################
+############################################## COMPILE SUMMA-Fortran ##############################################
+###################################################################################################################
+compile_noah:
+ $(FC) $(FLAGS_NOAH) -c $(NRUTIL) $(NOAHMP)
+
+# compile common routines
+compile_comm:
+ $(FC) $(FLAGS_COMM) -c $(COMM_ALL) $(INCLUDES)
+
+# compile SUMMA routines
+compile_summa:
+ $(FC) $(FLAGS_SUMMA) -c $(SUMMA_ALL) $(DRIVER) $(JOB_INTERFACE) $(FILEACCESS_INTERFACE) $(HRU_INTERFACE) $(GRU_INTERFACE) $(INCLUDES)
+
+# generate library
+link:
+ $(FC) -shared *.o -o libsumma.so
+ mv libsumma.so $(ROOT_DIR)/bin
+
+# Remove object files
+clean_fortran:
+ rm -f *.o *.mod soil_veg_gen_parm__genmod.f90
+###################################################################################################################
+############################################## COMPILE SUMMA-Fortran ##############################################
+###################################################################################################################
+
+
+###################################################################################################################
+################################################ COMPILE SUMMA-C++ ################################################
+###################################################################################################################
+compile_globals:
+ $(CC) $(FLAGS_ACTORS) -c $(GLOBAL) $(TIMEINFO) $(AUXILARY) $(SETTINGS_FILES) $(GLOBAL_INCLUDES) $(SUMMA_ACTOR_INCLUDES)
+
+compile_gru_actor:
+ $(CC) $(FLAGS_ACTORS) -c $(GRU_ACTOR) $(HRU_ACTOR_INCLUDES) $(GRU_ACTOR_INCLUDES) $(GLOBAL_INCLUDES) $(ACTORS_INCLUDES) $(SUMMA_ACTOR_INCLUDES)
+
+compile_hru_actor:
+ $(CC) $(FLAGS_ACTORS) -c $(HRU_ACTOR) $(HRU_ACTOR_INCLUDES) $(GLOBAL_INCLUDES) $(ACTORS_INCLUDES) $(SUMMA_ACTOR_INCLUDES)
+
+compile_file_access_actor:
+ $(CC) $(FLAGS_ACTORS) -c $(FILE_ACCESS_ACTOR) $(FORCING_FILE_INFO) $(OUTPUT_MANAGER) \
+ $(FILE_ACCESS_ACTOR_INCLUDES) $(GLOBAL_INCLUDES) $(ACTORS_INCLUDES) $(SUMMA_ACTOR_INCLUDES)
+
+compile_job_actor:
+ $(CC) $(FLAGS_ACTORS) -c $(JOB_ACTOR) $(GRUinfo) $(JOB_ACTOR_INCLUDES) $(GLOBAL_INCLUDES) $(ACTORS_INCLUDES) \
+ $(FILE_ACCESS_ACTOR_INCLUDES) $(HRU_ACTOR_INCLUDES) $(GRU_ACTOR_INCLUDES) $(SUMMA_ACTOR_INCLUDES)
+
+compile_summa_actor:
+ $(CC) $(FLAGS_ACTORS) -c $(SUMMA_ACTOR) $(SUMMA_ACTOR_INCLUDES) $(GLOBAL_INCLUDES) $(ACTORS_INCLUDES) \
+ $(JOB_ACTOR_INCLUDES) $(SUMMA_ACTOR_INCLUDES)
+
+compile_summa_client:
+ $(CC) $(FLAGS_ACTORS) -c $(SUMMA_CLIENT) $(SUMMA_ACTOR_INCLUDES) $(GLOBAL_INCLUDES)
+
+compile_summa_server:
+ $(CC) $(FLAGS_ACTORS) -c $(SUMMA_SERVER) $(SUMMA_BACKUP_SERVER) $(BATCH) $(CLIENT) $(BATCH_CONTAINER) $(CLIENT_CONTAINER) \
+ $(SUMMA_ACTOR_INCLUDES) $(GLOBAL_INCLUDES)
+
+compile_main:
+ $(CC) $(FLAGS_ACTORS) -c $(MAIN) $(GLOBAL_INCLUDES) $(SUMMA_ACTOR_INCLUDES) $(JOB_ACTOR_INCLUDES)
+
+link_cpp:
+ $(CC) $(FLAGS_ACTORS) -Wl,-rpath='/Summa-Actors/bin:/usr/local/lib:/code/sundials/instdir/lib' -o summaMain *.o $(ACTORS_LIBRARIES)
+ mv summaMain $(ROOT_DIR)/bin
+
+clean_cpp:
+ rm *.o
+###################################################################################################################
+################################################ COMPILE SUMMA-C++ ################################################
+###################################################################################################################
+
+
+###################################################################################################################
+################################################## COMPILE TESTS ##################################################
+###################################################################################################################
+actors_test:
+ $(CC) $(FLAGS_ACTORS) -c $(ACTOR_TEST) -std=c++17 $(ACTORS_INCLUDES)
+
+actors_testLink:
+ $(CC) -o summaTest *.o $(ACTORS_LIBRARIES)
+
+clean_lib:
+ rm *.so
+
+
+
+
diff --git a/build/makefile_sundials b/build/makefile_sundials
index 84e5021..17003f7 100644
--- a/build/makefile_sundials
+++ b/build/makefile_sundials
@@ -147,7 +147,7 @@ SUMMA_SOLVER= \
varSubstep.f90 \
sundials/varSubstepSundials.f90 \
opSplittin.f90 \
- coupled_em.f90
+ sundials/coupled_em.f90
SOLVER = $(patsubst %, $(ENGINE_DIR)/%, $(SUMMA_SOLVER))
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 b37b174..1a1cad0 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
@@ -21,7 +21,7 @@ module cppwrap_fileAccess
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
- USE summaActors_FileManager,only:LOCALPARAM_INFO,BASINPARAM_INFO ! files defining the default values and constraints for model parameters
+ USE summaFileManager,only:LOCALPARAM_INFO,BASINPARAM_INFO ! files defining the default values and constraints for model parameters
USE globalData,only:mpar_meta,bpar_meta ! parameter metadata structures
USE read_pinit_module,only:read_pinit ! module to read initial model parameter values
@@ -49,8 +49,8 @@ end subroutine read_pinit_C
subroutine read_vegitationTables(err) bind(C, name="read_vegitationTables")
USE SummaActors_setup,only:SOIL_VEG_GEN_PARM
USE module_sf_noahmplsm,only:read_mp_veg_parameters ! module to read NOAH vegetation tables
- USE summaActors_FileManager,only:SETTINGS_PATH ! define path to settings files (e.g., parameters, soil and veg. tables)
- USE summaActors_FileManager,only:GENPARM,VEGPARM,SOILPARM,MPTABLE ! files defining the noah tables
+ USE summaFileManager,only:SETTINGS_PATH ! define path to settings files (e.g., parameters, soil and veg. tables)
+ USE summaFileManager,only:GENPARM,VEGPARM,SOILPARM,MPTABLE ! files defining the noah tables
USE globalData,only:model_decisions ! model decision structure
USE var_lookup,only:iLookDECISIONS ! look-up values for model decisions
diff --git a/build/source/actors/file_access_actor/fortran_code/read_attribute.f90 b/build/source/actors/file_access_actor/fortran_code/read_attribute.f90
index 878c811..ff82536 100644
--- a/build/source/actors/file_access_actor/fortran_code/read_attribute.f90
+++ b/build/source/actors/file_access_actor/fortran_code/read_attribute.f90
@@ -55,8 +55,8 @@ subroutine openAttributeFile(attr_ncid, err) bind(C, name="openAttributeFile")
USE netcdf
USE netcdf_util_module,only:nc_file_open ! open netcdf file
! Attribute File
- USE summaActors_FileManager,only:SETTINGS_PATH ! define path to settings files (e.g., parameters, soil and veg. tables)
- USE summaActors_FileManager,only:LOCAL_ATTRIBUTES ! name of model initial attributes file
+ USE summaFileManager,only:SETTINGS_PATH ! define path to settings files (e.g., parameters, soil and veg. tables)
+ USE summaFileManager,only:LOCAL_ATTRIBUTES ! name of model initial attributes file
implicit none
integer(c_int),intent(out) :: attr_ncid
integer(c_int),intent(out) :: err
@@ -130,8 +130,8 @@ subroutine readAttributeFromNetCDF(ncid, index_gru, index_hru, num_var, &
USE globalData,only:attr_meta,type_meta,id_meta ! metadata structures
USE get_ixname_module,only:get_ixAttr,get_ixType,get_ixId ! access function to find index of elements in structure
! Information to make up the attributes file
- USE summaActors_FileManager,only:SETTINGS_PATH ! define path to settings files (e.g., parameters, soil and veg. tables)
- USE summaActors_FileManager,only:LOCAL_ATTRIBUTES ! name of model initial attributes file
+ USE summaFileManager,only:SETTINGS_PATH ! define path to settings files (e.g., parameters, soil and veg. tables)
+ USE summaFileManager,only:LOCAL_ATTRIBUTES ! name of model initial attributes file
! Fortran Data Type Structures
USE data_types,only:var_d, var_i, var_i8
implicit none
diff --git a/build/source/actors/file_access_actor/fortran_code/read_forcing.f90 b/build/source/actors/file_access_actor/fortran_code/read_forcing.f90
index be14d4b..18d639e 100644
--- a/build/source/actors/file_access_actor/fortran_code/read_forcing.f90
+++ b/build/source/actors/file_access_actor/fortran_code/read_forcing.f90
@@ -20,7 +20,7 @@ USE globalData,only:integerMissing ! integer missing value
USE var_lookup,only:iLookTIME,iLookFORCE ! named variables to define structure elements
-USE summaActors_FileManager,only:FORCING_PATH ! path of the forcing data file
+USE summaFileManager,only:FORCING_PATH ! path of the forcing data file
USE netcdf_util_module,only:nc_file_close ! close netcdf file
@@ -161,7 +161,7 @@ subroutine openForcingFile(forcFileInfo,iFile,infile,ncId,err,message)
USE globalData,only:utcTime ! all times in UTC (timeOffset = longitude/15. hours)
USE globalData,only:localTime ! all times local (timeOffset = 0)
USE globalData,only:refJulday_data
- USE summaActors_filemanager,only:NC_TIME_ZONE
+ USE summafilemanager,only:NC_TIME_ZONE
! dummy variables
type(file_info),intent(inout) :: forcFileInfo(:)
integer(i4b),intent(in) :: iFile ! index of current forcing file in forcing file list
diff --git a/build/source/actors/file_access_actor/fortran_code/read_param.f90 b/build/source/actors/file_access_actor/fortran_code/read_param.f90
index 52263b1..9fa29f3 100644
--- a/build/source/actors/file_access_actor/fortran_code/read_param.f90
+++ b/build/source/actors/file_access_actor/fortran_code/read_param.f90
@@ -59,8 +59,8 @@ subroutine openParamFile(param_ncid, param_file_exists, err) bind(C, name="openP
USE netcdf
USE netcdf_util_module,only:nc_file_open
! Parameters File
- USE summaActors_FileManager,only:SETTINGS_PATH ! path for metadata files
- USE summaActors_FileManager,only:PARAMETER_TRIAL ! file with parameter trial values
+ USE summaFileManager,only:SETTINGS_PATH ! path for metadata files
+ USE summaFileManager,only:PARAMETER_TRIAL ! file with parameter trial values
implicit none
integer(c_int),intent(out) :: param_ncid
diff --git a/build/source/actors/hru_actor/fortran_code/hru_actor.f90 b/build/source/actors/hru_actor/fortran_code/hru_actor.f90
index 28ecae6..1f14acb 100644
--- a/build/source/actors/hru_actor/fortran_code/hru_actor.f90
+++ b/build/source/actors/hru_actor/fortran_code/hru_actor.f90
@@ -70,7 +70,7 @@ subroutine setTimeZoneOffset(iFile, tmZoneOffsetFracDay, err) bind(C, name="setT
USE time_utils_module,only:extractTime ! extract time info from units string
USE time_utils_module,only:fracDay ! compute fractional day
- USE summaActors_filemanager,only:NC_TIME_ZONE
+ USE summafilemanager,only:NC_TIME_ZONE
implicit none
integer(c_int),intent(in) :: iFile
diff --git a/build/source/driver/summaActors_restart.f90 b/build/source/driver/summaActors_restart.f90
index 16ee61b..e20a2a3 100755
--- a/build/source/driver/summaActors_restart.f90
+++ b/build/source/driver/summaActors_restart.f90
@@ -77,9 +77,9 @@ contains
! global data structures
USE globalData,only:model_decisions ! model decision structure
! file paths
- USE summaActors_FileManager,only:SETTINGS_PATH ! path to settings files (e.g., Noah vegetation tables)
- USE summaActors_FileManager,only:STATE_PATH ! optional path to state/init. condition files (defaults to SETTINGS_PATH)
- USE summaActors_FileManager,only:MODEL_INITCOND ! name of model initial conditions file
+ USE summaFileManager,only:SETTINGS_PATH ! path to settings files (e.g., Noah vegetation tables)
+ USE summaFileManager,only:STATE_PATH ! optional path to state/init. condition files (defaults to SETTINGS_PATH)
+ USE summaFileManager,only:MODEL_INITCOND ! name of model initial conditions file
! timing variables
USE globalData,only:startRestart,endRestart ! date/time for the start and end of reading model restart files
USE globalData,only:elapsedRestart ! elapsed time to read model restart files
diff --git a/build/source/dshare/popMetadat.f90 b/build/source/dshare/popMetadat.f90
index 14dc59b..814bcd1 100755
--- a/build/source/dshare/popMetadat.f90
+++ b/build/source/dshare/popMetadat.f90
@@ -737,8 +737,8 @@ end subroutine popMetadat
subroutine read_output_file(err,message)
USE netcdf
! to get name of output control file from user
- USE summaActors_FileManager,only:SETTINGS_PATH ! path for metadata files
- USE summaActors_FileManager,only:OUTPUT_CONTROL ! file with output controls
+ USE summaFileManager,only:SETTINGS_PATH ! path for metadata files
+ USE summaFileManager,only:OUTPUT_CONTROL ! file with output controls
! some dimensional parameters
USE globalData, only:outFreq ! output frequencies
diff --git a/build/source/engine/coupled_em.f90 b/build/source/engine/coupled_em.f90
old mode 100755
new mode 100644
index da96f26..a2ed723
--- a/build/source/engine/coupled_em.f90
+++ b/build/source/engine/coupled_em.f90
@@ -65,8 +65,8 @@ USE globalData,only:globalPrintFlag ! the global print flag
! look-up values for the numerical method
USE mDecisions_module,only: &
- bEuler, & ! home-grown backward Euler solution with long time steps
- sundials ! SUNDIALS/IDA solution
+ bEuler, & ! home-grown backward Euler solution with long time steps
+ sundials ! SUNDIALS/IDA solution
! look-up values for the maximum interception capacity
USE mDecisions_module,only: &
@@ -96,31 +96,31 @@ real(dp),parameter :: dx=1.e-6_dp ! finite difference increment
contains
- ! ************************************************************************************************
- ! public subroutine coupled_em: run the coupled energy-mass model for one timestep
- ! ************************************************************************************************
+ ! ************************************************************************************************
+ ! public subroutine coupled_em: run the coupled energy-mass model for one timestep
+ ! ************************************************************************************************
subroutine coupled_em(&
- ! model control
- indxHRU, & ! intent(in): hruId
- dt_init, & ! intent(inout): used to initialize the size of the sub-step
- dt_init_factor, & ! Used to adjust the length of the timestep in the event of a failure
- computeVegFlux, & ! intent(inout): flag to indicate if we are computing fluxes over vegetation (.false. means veg is buried with snow)
- ! data structures (input)
- type_data, & ! intent(in): local classification of soil veg etc. for each HRU
- attr_data, & ! intent(in): local attributes for each HRU
- forc_data, & ! intent(in): model forcing data
- mpar_data, & ! intent(in): model parameters
- bvar_data, & ! intent(in): basin-average variables
- lookup_data, & ! intent(in): lookup tables
- ! data structures (input-output)
- indx_data, & ! intent(inout): model indices
- prog_data, & ! intent(inout): prognostic variables for a local HRU
- diag_data, & ! intent(inout): diagnostic variables for a local HRU
- flux_data, & ! intent(inout): model fluxes for a local HRU
- fracJulDay, &
- yearLength, &
- ! error control
- err,message) ! intent(out): error control
+ ! model control
+ indxHRU, & ! intent(in): hruId
+ dt_init, & ! intent(inout): used to initialize the size of the sub-step
+ dt_init_factor, & ! Used to adjust the length of the timestep in the event of a failure
+ computeVegFlux, & ! intent(inout): flag to indicate if we are computing fluxes over vegetation (.false. means veg is buried with snow)
+ ! data structures (input)
+ type_data, & ! intent(in): local classification of soil veg etc. for each HRU
+ attr_data, & ! intent(in): local attributes for each HRU
+ forc_data, & ! intent(in): model forcing data
+ mpar_data, & ! intent(in): model parameters
+ bvar_data, & ! intent(in): basin-average variables
+ lookup_data, & ! intent(in): lookup tables
+ ! data structures (input-output)
+ indx_data, & ! intent(inout): model indices
+ prog_data, & ! intent(inout): prognostic variables for a local HRU
+ diag_data, & ! intent(inout): diagnostic variables for a local HRU
+ flux_data, & ! intent(inout): model fluxes for a local HRU
+ fracJulDay, &
+ yearLength, &
+ ! error control
+ err,message) ! intent(out): error control
! structure allocations
USE allocspace_module,only:allocLocal ! allocate local data structures
USE allocspace_module,only:resizeData ! clone a data structure
@@ -250,7 +250,7 @@ subroutine coupled_em(&
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)
@@ -289,7 +289,7 @@ subroutine coupled_em(&
nLayers = nSnow + nSoil
- ! create temporary data structures for prognostic variables
+ ! 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
@@ -333,7 +333,7 @@ subroutine coupled_em(&
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(&
@@ -436,7 +436,7 @@ subroutine coupled_em(&
print*, message
return
end if
-
+
! check
if(computeVegFlux)then
@@ -662,7 +662,7 @@ subroutine coupled_em(&
! save/recover copies of index variables
do iVar=1,size(indx_data%var)
- !print*, 'indx_meta(iVar)%varname = ', trim(indx_meta(iVar)%varname)
+ !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(:)
@@ -805,16 +805,16 @@ subroutine coupled_em(&
! (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)
- prog_data%var(iLookPROG%scalarSfcMeltPond)%dat(1), & ! intent(inout): surface melt pond (kg m-2)
- ! input/output: properties of the upper-most soil layer
- prog_data%var(iLookPROG%mLayerTemp)%dat(nSnow+1), & ! intent(inout): surface layer temperature (K)
- prog_data%var(iLookPROG%mLayerDepth)%dat(nSnow+1), & ! intent(inout): surface layer depth (m)
- 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
+ ! 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)
+ prog_data%var(iLookPROG%scalarSfcMeltPond)%dat(1), & ! intent(inout): surface melt pond (kg m-2)
+ ! input/output: properties of the upper-most soil layer
+ prog_data%var(iLookPROG%mLayerTemp)%dat(nSnow+1), & ! intent(inout): surface layer temperature (K)
+ prog_data%var(iLookPROG%mLayerDepth)%dat(nSnow+1), & ! intent(inout): surface layer depth (m)
+ 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)
@@ -850,7 +850,6 @@ subroutine coupled_em(&
diag_data, & ! intent(inout): model diagnostic variables for a local HRU
flux_data, & ! intent(inout): model fluxes for a local HRU
bvar_data, & ! intent(in): model variables for the local basin
- lookup_data, & ! intent(in): lookup tables
model_decisions, & ! intent(in): model decisions
! output: model control
dtMultiplier, & ! intent(out): substep multiplier (-)
@@ -1005,7 +1004,7 @@ subroutine coupled_em(&
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))
@@ -1013,16 +1012,16 @@ subroutine coupled_em(&
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 *****************
@@ -1210,11 +1209,11 @@ subroutine coupled_em(&
! 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)
+ ! 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)
@@ -1227,17 +1226,17 @@ subroutine coupled_em(&
flux_data%var(averageFlux_meta(iVar)%ixParent)%dat(:) = flux_mean%var(iVar)%dat(:)
end do
- ! ***********************************************************************************************************************************
- ! ***********************************************************************************************************************************
- ! ***********************************************************************************************************************************
- ! ***********************************************************************************************************************************
+ ! ***********************************************************************************************************************************
+ ! ***********************************************************************************************************************************
+ ! ***********************************************************************************************************************************
+ ! ***********************************************************************************************************************************
- ! ---
- ! *** balance checks...
- ! ---------------------
+ ! ---
+ ! *** balance checks...
+ ! ---------------------
- ! save the average compression and melt pond storage in the data structures
- prog_data%var(iLookPROG%scalarSfcMeltPond)%dat(1) = sfcMeltPond
+ ! 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(&
@@ -1450,67 +1449,67 @@ subroutine coupled_em(&
end subroutine coupled_em
- ! *********************************************************************************************************
- ! private subroutine implctMelt: compute melt of the "snow without a layer"
- ! *********************************************************************************************************
- subroutine implctMelt(&
- ! input/output: integrated snowpack properties
- scalarSWE, & ! intent(inout): snow water equivalent (kg m-2)
- scalarSnowDepth, & ! intent(inout): snow depth (m)
- scalarSfcMeltPond, & ! intent(inout): surface melt pond (kg m-2)
- ! input/output: properties of the upper-most soil layer
- soilTemp, & ! intent(inout): surface layer temperature (K)
- soilDepth, & ! intent(inout): surface layer depth (m)
- soilHeatcap, & ! intent(inout): surface layer volumetric heat capacity (J m-3 K-1)
- ! output: error control
- err,message ) ! intent(out): error control
- implicit none
- ! input/output: integrated snowpack properties
- real(dp),intent(inout) :: scalarSWE ! snow water equivalent (kg m-2)
- real(dp),intent(inout) :: scalarSnowDepth ! snow depth (m)
- real(dp),intent(inout) :: scalarSfcMeltPond ! surface melt pond (kg m-2)
- ! input/output: properties of the upper-most soil layer
- real(dp),intent(inout) :: soilTemp ! surface layer temperature (K)
- real(dp),intent(inout) :: soilDepth ! surface layer depth (m)
- real(dp),intent(inout) :: soilHeatcap ! surface layer volumetric heat capacity (J m-3 K-1)
- ! output: error control
- integer(i4b),intent(out) :: err ! error code
- character(*),intent(out) :: message ! error message
- ! local variables
- real(dp) :: nrgRequired ! energy required to melt all the snow (J m-2)
- real(dp) :: nrgAvailable ! energy available to melt the snow (J m-2)
- real(dp) :: snwDensity ! snow density (kg m-3)
- ! initialize error control
- err=0; message='implctMelt/'
-
- if(scalarSWE > 0._dp)then
+ ! *********************************************************************************************************
+ ! private subroutine implctMelt: compute melt of the "snow without a layer"
+ ! *********************************************************************************************************
+ subroutine implctMelt(&
+ ! input/output: integrated snowpack properties
+ scalarSWE, & ! intent(inout): snow water equivalent (kg m-2)
+ scalarSnowDepth, & ! intent(inout): snow depth (m)
+ scalarSfcMeltPond, & ! intent(inout): surface melt pond (kg m-2)
+ ! input/output: properties of the upper-most soil layer
+ soilTemp, & ! intent(inout): surface layer temperature (K)
+ soilDepth, & ! intent(inout): surface layer depth (m)
+ soilHeatcap, & ! intent(inout): surface layer volumetric heat capacity (J m-3 K-1)
+ ! output: error control
+ err,message ) ! intent(out): error control
+ implicit none
+ ! input/output: integrated snowpack properties
+ real(dp),intent(inout) :: scalarSWE ! snow water equivalent (kg m-2)
+ real(dp),intent(inout) :: scalarSnowDepth ! snow depth (m)
+ real(dp),intent(inout) :: scalarSfcMeltPond ! surface melt pond (kg m-2)
+ ! input/output: properties of the upper-most soil layer
+ real(dp),intent(inout) :: soilTemp ! surface layer temperature (K)
+ real(dp),intent(inout) :: soilDepth ! surface layer depth (m)
+ real(dp),intent(inout) :: soilHeatcap ! surface layer volumetric heat capacity (J m-3 K-1)
+ ! output: error control
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! local variables
+ real(dp) :: nrgRequired ! energy required to melt all the snow (J m-2)
+ real(dp) :: nrgAvailable ! energy available to melt the snow (J m-2)
+ real(dp) :: snwDensity ! snow density (kg m-3)
+ ! initialize error control
+ err=0; message='implctMelt/'
+
+ if(scalarSWE > 0._dp)then
! only melt if temperature of the top soil layer is greater than Tfreeze
if(soilTemp > Tfreeze)then
- ! compute the energy required to melt all the snow (J m-2)
- nrgRequired = scalarSWE*LH_fus
- ! compute the energy available to melt the snow (J m-2)
- nrgAvailable = soilHeatcap*(soilTemp - Tfreeze)*soilDepth
- ! compute the snow density (not saved)
- snwDensity = scalarSWE/scalarSnowDepth
- ! compute the amount of melt, and update SWE (kg m-2)
- if(nrgAvailable > nrgRequired)then
+ ! compute the energy required to melt all the snow (J m-2)
+ nrgRequired = scalarSWE*LH_fus
+ ! compute the energy available to melt the snow (J m-2)
+ nrgAvailable = soilHeatcap*(soilTemp - Tfreeze)*soilDepth
+ ! compute the snow density (not saved)
+ snwDensity = scalarSWE/scalarSnowDepth
+ ! compute the amount of melt, and update SWE (kg m-2)
+ if(nrgAvailable > nrgRequired)then
scalarSfcMeltPond = scalarSWE
scalarSWE = 0._dp
- else
+ else
scalarSfcMeltPond = nrgAvailable/LH_fus
scalarSWE = scalarSWE - scalarSfcMeltPond
- end if
- ! update depth
- scalarSnowDepth = scalarSWE/snwDensity
- ! update temperature of the top soil layer (K)
- soilTemp = soilTemp - (LH_fus*scalarSfcMeltPond/soilDepth)/soilHeatcap
+ end if
+ ! update depth
+ scalarSnowDepth = scalarSWE/snwDensity
+ ! update temperature of the top soil layer (K)
+ soilTemp = soilTemp - (LH_fus*scalarSfcMeltPond/soilDepth)/soilHeatcap
else ! melt is zero if the temperature of the top soil layer is less than Tfreeze
- scalarSfcMeltPond = 0._dp ! kg m-2
+ scalarSfcMeltPond = 0._dp ! kg m-2
end if ! (if the temperature of the top soil layer is greater than Tfreeze)
- else ! melt is zero if the "snow without a layer" does not exist
+ else ! melt is zero if the "snow without a layer" does not exist
scalarSfcMeltPond = 0._dp ! kg m-2
- end if ! (if the "snow without a layer" exists)
+ end if ! (if the "snow without a layer" exists)
- end subroutine implctMelt
+ end subroutine implctMelt
end module coupled_em_module
diff --git a/build/source/engine/derivforce.f90 b/build/source/engine/derivforce.f90
index a2d75e7..d4909d6 100755
--- a/build/source/engine/derivforce.f90
+++ b/build/source/engine/derivforce.f90
@@ -69,7 +69,7 @@ contains
USE conv_funcs_module,only:SPHM2RELHM,RELHM2SPHM,WETBULBTMP ! conversion functions
USE snow_utils_module,only:fracliquid,templiquid ! functions to compute temperature/liquid water
USE time_utils_module,only:compcalday ! convert julian day to calendar date
- USE summaActors_FileManager,only: NC_TIME_ZONE ! time zone option from control file
+ USE summaFileManager,only: NC_TIME_ZONE ! time zone option from control file
! compute derived forcing data variables
implicit none
! input variables
diff --git a/build/source/engine/ffile_info.f90 b/build/source/engine/ffile_info.f90
index 863c545..02f598c 100755
--- a/build/source/engine/ffile_info.f90
+++ b/build/source/engine/ffile_info.f90
@@ -44,9 +44,9 @@ subroutine ffile_info(indxGRU,handle_forcFileInfo,num_forcing_files,err) bind(C,
USE ascii_util_module,only:linewidth
USE netcdf_util_module,only:nc_file_open ! open netCDF file
USE netcdf_util_module,only:netcdf_err ! netcdf error handling function
- USE summaActors_FileManager,only:SETTINGS_PATH ! path for metadata files
- USE summaActors_FileManager,only:FORCING_PATH ! path for forcing files
- USE summaActors_FileManager,only:FORCING_FILELIST ! list of model forcing files
+ USE summaFileManager,only:SETTINGS_PATH ! path for metadata files
+ USE summaFileManager,only:FORCING_PATH ! path for forcing files
+ USE summaFileManager,only:FORCING_FILELIST ! list of model forcing files
USE globalData,only:forc_meta ! forcing metadata
USE get_ixname_module,only:get_ixtime,get_ixforce ! identify index of named variable
USE ascii_util_module,only:get_vlines ! get a vector of non-comment lines
@@ -55,7 +55,7 @@ subroutine ffile_info(indxGRU,handle_forcFileInfo,num_forcing_files,err) bind(C,
USE globalData,only:time_meta
USE allocspace_module,only:allocLocal
USE time_utils_module,only:extractTime ! extract time info from units string
- USE summaActors_FileManager,only: SIM_START_TM, SIM_END_TM! time info from control file module
+ USE summaFileManager,only: SIM_START_TM, SIM_END_TM! time info from control file module
USE var_lookup,only:iLookTIME ! named variables that identify indices in the time structures
diff --git a/build/source/engine/mDecisions.f90 b/build/source/engine/mDecisions.f90
index 595f37c..2a850ad 100755
--- a/build/source/engine/mDecisions.f90
+++ b/build/source/engine/mDecisions.f90
@@ -182,7 +182,7 @@ subroutine mDecisions(num_steps,err) bind(C, name='mDecisions')
USE time_utils_module,only:extractTime ! extract time info from units string
USE time_utils_module,only:compjulday ! compute the julian day
USE time_utils_module,only:fracDay ! compute fractional day
- USE summaActors_FileManager,only: SIM_START_TM, SIM_END_TM ! time info from control file module
+ USE summaFileManager,only: SIM_START_TM, SIM_END_TM ! time info from control file module
implicit none
! define output
@@ -677,8 +677,8 @@ subroutine readoption(err,message)
USE ascii_util_module,only:file_open ! open file
USE ascii_util_module,only:linewidth ! max character number for one line
USE ascii_util_module,only:get_vlines ! get a vector of non-comment lines
- USE summaActors_FileManager,only:SETTINGS_PATH ! path for metadata files
- USE summaActors_FileManager,only:M_DECISIONS ! definition of modeling options
+ USE summaFileManager,only:SETTINGS_PATH ! path for metadata files
+ USE summaFileManager,only:M_DECISIONS ! definition of modeling options
USE get_ixname_module,only:get_ixdecisions ! identify index of named variable
USE globalData,only:model_decisions ! model decision structure
implicit none
diff --git a/build/source/engine/read_dimension.f90 b/build/source/engine/read_dimension.f90
index f8da759..404c1b1 100644
--- a/build/source/engine/read_dimension.f90
+++ b/build/source/engine/read_dimension.f90
@@ -40,8 +40,8 @@ subroutine read_dimension(numGRUs,numHRUs,startGRU,err) bind(C, name="readDimens
USE globalData,only:gru_struc ! gru->hru mapping structure
USE globalData,only:index_map ! hru->gru mapping structure
! file paths for attribute file
- USE summaActors_FileManager,only:SETTINGS_PATH ! define path to settings files (e.g., parameters, soil and veg. tables)
- USE summaActors_FileManager,only:LOCAL_ATTRIBUTES ! name of model initial attributes file
+ USE summaFileManager,only:SETTINGS_PATH ! define path to settings files (e.g., parameters, soil and veg. tables)
+ USE summaFileManager,only:LOCAL_ATTRIBUTES ! name of model initial attributes file
implicit none
diff --git a/build/source/engine/read_pinit.f90 b/build/source/engine/read_pinit.f90
index 5c36ff4..2a0b350 100755
--- a/build/source/engine/read_pinit.f90
+++ b/build/source/engine/read_pinit.f90
@@ -36,7 +36,7 @@ contains
! ************************************************************************************************
subroutine read_pinit(filenm,isLocal,mpar_meta,parFallback,err,message)
! used to read metadata on the forcing data file
- USE summaActors_FileManager,only:SETTINGS_PATH ! path for input parameter and other configuration files
+ USE summaFileManager,only:SETTINGS_PATH ! path for input parameter and other configuration files
USE ascii_util_module,only:file_open ! open ascii file
USE ascii_util_module,only:split_line ! extract the list of variable names from the character string
USE data_types,only:var_info ! data type for metadata
diff --git a/build/source/engine/sundials/coupled_em.f90 b/build/source/engine/sundials/coupled_em.f90
new file mode 100755
index 0000000..da96f26
--- /dev/null
+++ b/build/source/engine/sundials/coupled_em.f90
@@ -0,0 +1,1516 @@
+! SUMMA - Structure for Unifying Multiple Modeling Alternatives
+! Copyright (C) 2014-2020 NCAR/RAL; University of Saskatchewan; University of Washington
+!
+! This file is part of SUMMA
+!
+! For more information see: http://www.ral.ucar.edu/projects/summa
+!
+! This program is free software: you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License
+! along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+module coupled_em_module
+
+! numerical recipes data types
+USE nrtype
+
+! physical constants
+USE multiconst,only:&
+ Tfreeze, & ! temperature at freezing (K)
+ LH_fus, & ! latent heat of fusion (J kg-1)
+ LH_sub, & ! latent heat of sublimation (J kg-1)
+ iden_ice, & ! intrinsic density of ice (kg m-3)
+ iden_water ! intrinsic density of liquid water (kg m-3)
+
+! data types
+USE data_types,only:&
+ var_i, & ! x%var(:) (i4b)
+ var_d, & ! x%var(:) (dp)
+ var_ilength, & ! x%var(:)%dat (i4b)
+ var_dlength, & ! x%var(:)%dat (dp)
+ zLookup
+
+! named variables for parent structures
+USE var_lookup,only:iLookDECISIONS ! named variables for elements of the decision structure
+USE var_lookup,only:iLookPROG ! named variables for structure elements
+USE var_lookup,only:iLookDIAG ! named variables for structure elements
+USE var_lookup,only:iLookFLUX ! named variables for structure elements
+USE var_lookup,only:iLookPARAM ! named variables for structure elements
+USE var_lookup,only:iLookINDEX ! named variables for structure elements
+USE globalData,only:iname_snow ! named variables for snow
+USE globalData,only:iname_soil ! named variables for soil
+
+! named variables for child structures
+USE var_lookup,only:childFLUX_MEAN
+
+! metadata
+USE globalData,only:indx_meta ! metadata on the model index variables
+USE globalData,only:diag_meta ! metadata on the model diagnostic variables
+USE globalData,only:prog_meta ! metadata on the model prognostic variables
+USE globalData,only:averageFlux_meta ! metadata on the timestep-average model flux structure
+
+! global data
+USE globalData,only:data_step ! time step of forcing data (s)
+USE globalData,only:model_decisions ! model decision structure
+USE globalData,only:globalPrintFlag ! the global print flag
+
+! look-up values for the numerical method
+USE mDecisions_module,only: &
+ bEuler, & ! home-grown backward Euler solution with long time steps
+ sundials ! SUNDIALS/IDA solution
+
+! 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
+
+! look-up values for the groundwater parameterization
+USE mDecisions_module,only: &
+ qbaseTopmodel,& ! TOPMODEL-ish baseflow parameterization
+ bigBucket ,& ! a big bucket (lumped aquifer model)
+ noExplicit ! no explicit groundwater parameterization
+
+! look-up values for the spatial representation of groundwater
+USE mDecisions_module,only: &
+ localColumn ,& ! separate groundwater representation in each local soil column
+ singleBasin ! single groundwater store over the entire basin
+
+! privacy
+implicit none
+private
+public::coupled_em
+! algorithmic parameters
+real(dp),parameter :: valueMissing=-9999._dp ! missing value, used when diagnostic or state variables are undefined
+real(dp),parameter :: verySmall=1.e-6_dp ! used as an additive constant to check if substantial difference among real numbers
+real(dp),parameter :: mpe=1.e-6_dp ! prevents overflow error if division by zero
+real(dp),parameter :: dx=1.e-6_dp ! finite difference increment
+contains
+
+
+ ! ************************************************************************************************
+ ! public subroutine coupled_em: run the coupled energy-mass model for one timestep
+ ! ************************************************************************************************
+subroutine coupled_em(&
+ ! model control
+ indxHRU, & ! intent(in): hruId
+ dt_init, & ! intent(inout): used to initialize the size of the sub-step
+ dt_init_factor, & ! Used to adjust the length of the timestep in the event of a failure
+ computeVegFlux, & ! intent(inout): flag to indicate if we are computing fluxes over vegetation (.false. means veg is buried with snow)
+ ! data structures (input)
+ type_data, & ! intent(in): local classification of soil veg etc. for each HRU
+ attr_data, & ! intent(in): local attributes for each HRU
+ forc_data, & ! intent(in): model forcing data
+ mpar_data, & ! intent(in): model parameters
+ bvar_data, & ! intent(in): basin-average variables
+ lookup_data, & ! intent(in): lookup tables
+ ! data structures (input-output)
+ indx_data, & ! intent(inout): model indices
+ prog_data, & ! intent(inout): prognostic variables for a local HRU
+ diag_data, & ! intent(inout): diagnostic variables for a local HRU
+ flux_data, & ! intent(inout): model fluxes for a local HRU
+ fracJulDay, &
+ yearLength, &
+ ! error control
+ err,message) ! intent(out): error control
+ ! structure allocations
+ USE allocspace_module,only:allocLocal ! allocate local data structures
+ USE allocspace_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
+ USE computSnowDepth_module,only:computSnowDepth
+ ! look-up values for the numerical method
+ 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
+ type(zLookup),intent(in) :: lookup_data ! lookup tables
+ ! 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=.true.
+ 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/"
+
+ ! 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 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
+
+ ! 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)
+
+ ! 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)
+ 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)
+ 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)
+ 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)
+ 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 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'
+ print*,message
+ err=20;
+ return
+ 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
+
+ ! 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)
+
+ ! save SWE
+ oldSWE = prog_data%var(iLookPROG%scalarSWE)%dat(1)
+ !print*, 'nSnow = ', nSnow
+ !print*, 'oldSWE = ', oldSWE
+
+ ! *** 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))
+
+ ! 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)
+ 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.'
+ print*, message
+ err=20; return
+ 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'
+ 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
+ ! 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
+ (prog_data%var(iLookPROG%scalarCanopyTemp)%dat(1) < Tfreeze), & ! flag to denote if the canopy is frozen
+ varNotUsed1, & ! derivative in canopy liquid w.r.t. canopy temperature (kg m-2 K-1)
+ varNotUsed2, & ! fraction of liquid water on the canopy
+ prog_data%var(iLookPROG%scalarCanopyLiq)%dat(1), & ! canopy liquid water (kg m-2)
+ prog_data%var(iLookPROG%scalarCanopyIce)%dat(1), & ! canopy ice (kg m-2)
+ diag_data%var(iLookDIAG%scalarCanopyLiqMax)%dat(1), & ! maximum canopy liquid water (kg m-2)
+ diag_data%var(iLookDIAG%scalarCanopyLiqMax)%dat(1), & ! maximum canopy ice content (kg m-2)
+ mpar_data%var(iLookPARAM%canopyWettingFactor)%dat(1), & ! maximum wetted fraction of the canopy (-)
+ mpar_data%var(iLookPARAM%canopyWettingExp)%dat(1), & ! exponent in canopy wetting function (-)
+ ! output
+ diag_data%var(iLookDIAG%scalarCanopyWetFraction)%dat(1), & ! canopy wetted fraction (-)
+ 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)
+ 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
+
+
+ ! *** 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)
+ 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)
+ 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(&
+ ! input: derived parameters
+ canopyDepth, & ! intent(in): canopy depth (m)
+ ! input/output: data structures
+ mpar_data, & ! intent(in): model parameters
+ prog_data, & ! intent(inout): model prognostic variables for a local HRU
+ 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
+ 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)
+ 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
+
+
+ ! *** 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(&
+ ! input: control variables
+ computeVegFlux, & ! intent(in): flag to denote if computing the vegetation flux
+ canopyDepth, & ! intent(in): canopy depth (m)
+ ! input/output: data structures
+ mpar_data, & ! intent(in): model parameters
+ indx_data, & ! intent(in): model layer indices
+ prog_data, & ! intent(in): model prognostic variables for a local HRU
+ diag_data, & ! intent(inout): model diagnostic variables for a local HRU
+ ! 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(&
+ ! 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)
+ prog_data%var(iLookPROG%scalarSfcMeltPond)%dat(1), & ! intent(inout): surface melt pond (kg m-2)
+ ! input/output: properties of the upper-most soil layer
+ prog_data%var(iLookPROG%mLayerTemp)%dat(nSnow+1), & ! intent(inout): surface layer temperature (K)
+ prog_data%var(iLookPROG%mLayerDepth)%dat(nSnow+1), & ! intent(inout): surface layer depth (m)
+ 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)
+ print*, message
+ return
+ end if
+ end if ! nsnow == 0
+
+ ! *** solve model equations...
+ ! ----------------------------
+
+ ! save input step
+ dtSave = dt_sub
+ !write(*,'(a,1x,3(f12.5,1x))') trim(message)//'before opSplittin: dt_init, dt_sub, dt_solv = ', dt_init, dt_sub, dt_solv
+
+ ! get the new solution
+ call opSplittin(&
+ ! input: model control
+ nSnow, & ! intent(in): number of snow layers
+ nSoil, & ! intent(in): number of soil layers
+ nLayers, & ! intent(in): total number of layers
+ nState, & ! intent(in): total number of layers
+ dt_sub, & ! intent(in): length of the model sub-step
+ (nsub==1), & ! intent(in): logical flag to denote the first substep
+ computeVegFlux, & ! intent(in): logical flag to compute fluxes within the vegetation canopy
+ ! input/output: data structures
+ type_data, & ! intent(in): type of vegetation and soil
+ attr_data, & ! intent(in): spatial attributes
+ forc_data, & ! intent(in): model forcing data
+ mpar_data, & ! intent(in): model parameters
+ indx_data, & ! intent(inout): index data
+ 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
+ bvar_data, & ! intent(in): model variables for the local basin
+ lookup_data, & ! intent(in): lookup tables
+ model_decisions, & ! intent(in): model decisions
+ ! output: model control
+ dtMultiplier, & ! intent(out): substep multiplier (-)
+ tooMuchMelt, & ! intent(out): flag to denote that ice is insufficient to support melt
+ stepFailure, & ! intent(out): flag to denote that the coupled step failed
+ ixSolution, & ! intent(out): solution method used in this iteration
+ 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)
+ print*, message
+ return
+ end if
+
+
+ ! process the flag for too much melt
+ if(tooMuchMelt)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
+ ! 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
+
+ ! 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
+
+ ! update first step
+ firstSubStep=.false.
+
+ ! *** remove ice due to sublimation...
+ ! --------------------------------------------------------------
+ sublime: associate(&
+ scalarCanopySublimation => flux_data%var(iLookFLUX%scalarCanopySublimation)%dat(1), & ! sublimation from the vegetation canopy (kg m-2 s-1)
+ scalarSnowSublimation => flux_data%var(iLookFLUX%scalarSnowSublimation)%dat(1), & ! sublimation from the snow surface (kg m-2 s-1)
+ scalarLatHeatCanopyEvap => flux_data%var(iLookFLUX%scalarLatHeatCanopyEvap)%dat(1), & ! latent heat flux for evaporation from the canopy to the canopy air space (W m-2)
+ scalarSenHeatCanopy => flux_data%var(iLookFLUX%scalarSenHeatCanopy)%dat(1), & ! sensible heat flux from the canopy to the canopy air space (W m-2)
+ scalarLatHeatGround => flux_data%var(iLookFLUX%scalarLatHeatGround)%dat(1), & ! latent heat flux from ground surface below vegetation (W m-2)
+ scalarSenHeatGround => flux_data%var(iLookFLUX%scalarSenHeatGround)%dat(1), & ! sensible heat flux from ground surface below vegetation (W m-2)
+ scalarCanopyLiq => prog_data%var(iLookPROG%scalarCanopyLiq)%dat(1), & ! liquid water stored on the vegetation canopy (kg m-2)
+ scalarCanopyIce => prog_data%var(iLookPROG%scalarCanopyIce)%dat(1), & ! ice stored on the vegetation canopy (kg m-2)
+ 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
+
+ ! * compute change in canopy ice content due to sublimation...
+ ! ------------------------------------------------------------
+ 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
+
+ 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
+
+ ! 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)
+ 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)
+
+ end if ! sundials
+ end associate sublime
+
+
+ ! ****************************************************************************************************
+ ! *** END MAIN SOLVER ********************************************************************************
+ ! ****************************************************************************************************
+
+ ! increment sub-step
+ dt_solv = dt_solv + dt_sub
+
+ ! save the time step to initialize the subsequent step
+ if(dt_solv<data_step .or. nsub==1) dt_init = dt_sub
+
+ ! check
+ if(globalPrintFlag)&
+ write(*,'(a,1x,3(f18.5,1x))') 'dt_sub, dt_solv, data_step: ', dt_sub, dt_solv, data_step
+
+ ! check that we have completed the sub-step
+ if(dt_solv >= data_step-verySmall) then
+ exit substeps
+ endif
+
+ ! adjust length of the sub-step (make sure that we don't exceed the step)
+ dt_sub = data_step - dt_solv
+
+ end do substeps ! (sub-step loop)
+
+ ! *** add snowfall to the snowpack...
+ ! -----------------------------------
+
+ ! add new snowfall to the snowpack
+ ! NOTE: This needs to be done AFTER the call to canopySnow, since throughfall and unloading are computed in canopySnow
+ call newsnwfall(&
+ ! input: model control
+ data_step, & ! time step (seconds)
+ (nSnow > 0), & ! logical flag if snow layers exist
+ mpar_data%var(iLookPARAM%snowfrz_scale)%dat(1), & ! freeezing curve parameter for snow (K-1)
+ ! input: diagnostic scalar variables
+ diag_data%var(iLookDIAG%scalarSnowfallTemp)%dat(1), & ! computed temperature of fresh snow (K)
+ diag_data%var(iLookDIAG%scalarNewSnowDensity)%dat(1), & ! computed density of new snow (kg m-3)
+ flux_data%var(iLookFLUX%scalarThroughfallSnow)%dat(1), & ! throughfall of snow through the canopy (kg m-2 s-1)
+ flux_data%var(iLookFLUX%scalarCanopySnowUnloading)%dat(1), & ! unloading of snow from the canopy (kg m-2 s-1)
+ ! input/output: state variables
+ prog_data%var(iLookPROG%scalarSWE)%dat(1), & ! SWE (kg m-2)
+ prog_data%var(iLookPROG%scalarSnowDepth)%dat(1), & ! total snow depth (m)
+ prog_data%var(iLookPROG%mLayerTemp)%dat(1), & ! temperature of the top layer (K)
+ prog_data%var(iLookPROG%mLayerDepth)%dat(1), & ! depth of the top layer (m)
+ prog_data%var(iLookPROG%mLayerVolFracIce)%dat(1), & ! volumetric fraction of ice of the top layer (-)
+ prog_data%var(iLookPROG%mLayerVolFracLiq)%dat(1), & ! volumetric fraction of liquid water of the top layer (-)
+ ! output: error control
+ err,cmessage) ! error control
+ if(err/=0)then
+ err=30
+ message=trim(message)//trim(cmessage)
+ print*,message
+ return
+ end if
+
+ ! re-compute 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
+ !print*, 'SWE after snowfall = ', prog_data%var(iLookPROG%scalarSWE)%dat(1)
+
+ ! 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
+
+ ! 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
+
+ ! overwrite flux_data with flux_mean (returns timestep-average fluxes for scalar variables)
+ do iVar=1,size(averageFlux_meta)
+ 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...
+ ! ----------------------------
+
+ ! 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
+ print*, message
+ 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
+
+ ! end association of local variables with information in the data structures
+ end associate
+
+ ! end association to canopy depth
+ end associate canopy
+
+ ! Save the total soil water (Liquid+Ice)
+ diag_data%var(iLookDIAG%scalarTotalSoilWat)%dat(1) = balanceSoilWater1
+ ! save the surface temperature (just to make things easier to visualize)
+ prog_data%var(iLookPROG%scalarSurfaceTemp)%dat(1) = prog_data%var(iLookPROG%mLayerTemp)%dat(1)
+
+ ! overwrite flux data with the timestep-average value
+ if(.not.backwardsCompatibility)then
+ do iVar=1,size(flux_mean%var)
+ flux_data%var(averageFlux_meta(iVar)%ixParent)%dat = flux_mean%var(iVar)%dat
+ end do
+ 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
+
+
+ ! *********************************************************************************************************
+ ! private subroutine implctMelt: compute melt of the "snow without a layer"
+ ! *********************************************************************************************************
+ subroutine implctMelt(&
+ ! input/output: integrated snowpack properties
+ scalarSWE, & ! intent(inout): snow water equivalent (kg m-2)
+ scalarSnowDepth, & ! intent(inout): snow depth (m)
+ scalarSfcMeltPond, & ! intent(inout): surface melt pond (kg m-2)
+ ! input/output: properties of the upper-most soil layer
+ soilTemp, & ! intent(inout): surface layer temperature (K)
+ soilDepth, & ! intent(inout): surface layer depth (m)
+ soilHeatcap, & ! intent(inout): surface layer volumetric heat capacity (J m-3 K-1)
+ ! output: error control
+ err,message ) ! intent(out): error control
+ implicit none
+ ! input/output: integrated snowpack properties
+ real(dp),intent(inout) :: scalarSWE ! snow water equivalent (kg m-2)
+ real(dp),intent(inout) :: scalarSnowDepth ! snow depth (m)
+ real(dp),intent(inout) :: scalarSfcMeltPond ! surface melt pond (kg m-2)
+ ! input/output: properties of the upper-most soil layer
+ real(dp),intent(inout) :: soilTemp ! surface layer temperature (K)
+ real(dp),intent(inout) :: soilDepth ! surface layer depth (m)
+ real(dp),intent(inout) :: soilHeatcap ! surface layer volumetric heat capacity (J m-3 K-1)
+ ! output: error control
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! local variables
+ real(dp) :: nrgRequired ! energy required to melt all the snow (J m-2)
+ real(dp) :: nrgAvailable ! energy available to melt the snow (J m-2)
+ real(dp) :: snwDensity ! snow density (kg m-3)
+ ! initialize error control
+ err=0; message='implctMelt/'
+
+ if(scalarSWE > 0._dp)then
+ ! only melt if temperature of the top soil layer is greater than Tfreeze
+ if(soilTemp > Tfreeze)then
+ ! compute the energy required to melt all the snow (J m-2)
+ nrgRequired = scalarSWE*LH_fus
+ ! compute the energy available to melt the snow (J m-2)
+ nrgAvailable = soilHeatcap*(soilTemp - Tfreeze)*soilDepth
+ ! compute the snow density (not saved)
+ snwDensity = scalarSWE/scalarSnowDepth
+ ! compute the amount of melt, and update SWE (kg m-2)
+ if(nrgAvailable > nrgRequired)then
+ scalarSfcMeltPond = scalarSWE
+ scalarSWE = 0._dp
+ else
+ scalarSfcMeltPond = nrgAvailable/LH_fus
+ scalarSWE = scalarSWE - scalarSfcMeltPond
+ end if
+ ! update depth
+ scalarSnowDepth = scalarSWE/snwDensity
+ ! update temperature of the top soil layer (K)
+ soilTemp = soilTemp - (LH_fus*scalarSfcMeltPond/soilDepth)/soilHeatcap
+ else ! melt is zero if the temperature of the top soil layer is less than Tfreeze
+ scalarSfcMeltPond = 0._dp ! kg m-2
+ end if ! (if the temperature of the top soil layer is greater than Tfreeze)
+ else ! melt is zero if the "snow without a layer" does not exist
+ scalarSfcMeltPond = 0._dp ! kg m-2
+ end if ! (if the "snow without a layer" exists)
+
+ end subroutine implctMelt
+
+end module coupled_em_module
diff --git a/build/source/engine/vegPhenlgy.f90 b/build/source/engine/vegPhenlgy.f90
index 07e14c2..7850909 100755
--- a/build/source/engine/vegPhenlgy.f90
+++ b/build/source/engine/vegPhenlgy.f90
@@ -25,8 +25,6 @@ USE nrtype
! global variables
USE globalData,only:urbanVegCategory ! vegetation category for urban areas
-!USE globalData,only:fracJulday ! fractional julian days since the start of year
-! USE globalData,only:yearLength ! number of days in the current year
! provide access to the derived types to define the data structures
USE data_types,only:&
diff --git a/build/source/hookup/summaActors_FileManager.f90 b/build/source/hookup/summaActors_FileManager.f90
index 0039367..161e849 100755
--- a/build/source/hookup/summaActors_FileManager.f90
+++ b/build/source/hookup/summaActors_FileManager.f90
@@ -22,7 +22,7 @@
! Original version based on:
! (C) Copyright 2009-2010 --- Dmitri Kavetski and Martyn Clark --- All rights reserved
!******************************************************************
-MODULE summaActors_FileManager
+MODULE summaFileManager
USE, intrinsic :: iso_c_binding
use nrtype
implicit none
@@ -180,4 +180,4 @@ subroutine summa_SetTimesDirsAndFiles(file_manager,err) bind(C, name="setTimesDi
end subroutine summa_SetTimesDirsAndFiles
-END MODULE summaActors_FileManager
+END MODULE summaFileManager
diff --git a/build/source/netcdf/def_output.f90 b/build/source/netcdf/def_output.f90
index 26f8ca6..8ff8a25 100755
--- a/build/source/netcdf/def_output.f90
+++ b/build/source/netcdf/def_output.f90
@@ -89,7 +89,7 @@ subroutine def_output(handle_ncid,startGRU,nGRU,nHRU,err) bind(C, name='def_outp
! modules that are not globalData
USE var_lookup,only:maxVarFreq ! # of available output frequencies
USE get_ixname_module,only:get_freqName ! get name of frequency from frequency index
- USE summaActors_FileManager,only:OUTPUT_PATH,OUTPUT_PREFIX ! define output file
+ USE summaFileManager,only:OUTPUT_PATH,OUTPUT_PREFIX ! define output file
! ---------------------------------------------------------------------------------------
! * variables from C++
diff --git a/build/source/netcdf/read_icondActors.f90 b/build/source/netcdf/read_icondActors.f90
index 52bf5a1..42422e1 100644
--- a/build/source/netcdf/read_icondActors.f90
+++ b/build/source/netcdf/read_icondActors.f90
@@ -52,9 +52,9 @@ subroutine read_icond_nlayers(nGRU,err) bind(C, name="readIcondNLayers")
USE globalData,only:indx_meta
! file paths
- USE summaActors_FileManager,only:STATE_PATH ! optional path to state/init. condition files (defaults to SETTINGS_PATH)
- USE summaActors_FileManager,only:SETTINGS_PATH ! define path to settings files (e.g., parameters, soil and veg. tables)
- USE summaActors_FileManager,only:MODEL_INITCOND ! name of model initial conditions file
+ USE summaFileManager,only:STATE_PATH ! optional path to state/init. condition files (defaults to SETTINGS_PATH)
+ USE summaFileManager,only:SETTINGS_PATH ! define path to settings files (e.g., parameters, soil and veg. tables)
+ USE summaFileManager,only:MODEL_INITCOND ! name of model initial conditions file
implicit none
diff --git a/build/summa b/build/summa
new file mode 160000
index 0000000..fa9adf8
--- /dev/null
+++ b/build/summa
@@ -0,0 +1 @@
+Subproject commit fa9adf808229a45085defdc2bb8ef05836b9b3aa
--
GitLab