diff --git a/.github/workflows/linux.yml b/.github/workflows/linux.yml index ca455507d..deeba49c2 100644 --- a/.github/workflows/linux.yml +++ b/.github/workflows/linux.yml @@ -397,7 +397,7 @@ jobs: make -j make install fi - echo "RMM_FLAGS=-Drmm_ROOT=$PARFLOW_DEP_DIR/rmm" >> $GITHUB_ENV + echo "RMM_FLAGS=-DRMM_ROOT=$PARFLOW_DEP_DIR/rmm" >> $GITHUB_ENV - name: Umpire Install env: @@ -414,7 +414,7 @@ jobs: cmake --build build --parallel 4 cmake --install build fi - echo "UMPIRE_FLAGS=-Dumpire_ROOT=$PARFLOW_DEP_DIR/Umpire" >> $GITHUB_ENV + echo "UMPIRE_FLAGS=-DUMPIRE_ROOT=$PARFLOW_DEP_DIR/Umpire" >> $GITHUB_ENV - name: Kokkos Install env: diff --git a/.github/workflows/self_hosted.yml b/.github/workflows/self_hosted.yml index 9f5494e1b..329026ba8 100644 --- a/.github/workflows/self_hosted.yml +++ b/.github/workflows/self_hosted.yml @@ -48,7 +48,7 @@ jobs: cmake .. -DCMAKE_INSTALL_PREFIX=$GITHUB_WORKSPACE/depend/rmm -DBUILD_TESTS=OFF make -j make install - echo "RMM_FLAGS=-Drmm_ROOT=$GITHUB_WORKSPACE/depend/rmm" >> $GITHUB_ENV + echo "RMM_FLAGS=-DRMM_ROOT=$GITHUB_WORKSPACE/depend/rmm" >> $GITHUB_ENV - name: Install Umpire if: matrix.config.memory_manager == 'umpire' @@ -63,7 +63,7 @@ jobs: cmake .. -DCMAKE_INSTALL_PREFIX=$GITHUB_WORKSPACE/depend/Umpire -DENABLE_CUDA=On make -j make install - echo "UMPIRE_FLAGS=-Dumpire_ROOT=$GITHUB_WORKSPACE/depend/Umpire" >> $GITHUB_ENV + echo "UMPIRE_FLAGS=-DUMPIRE_ROOT=$GITHUB_WORKSPACE/depend/Umpire" >> $GITHUB_ENV - name: Configure ParFlow run: | diff --git a/CMakeLists.txt b/CMakeLists.txt index c2179b0ae..4f3c5b59a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -125,27 +125,27 @@ else(PARFLOW_ACCELERATOR_BACKEND STREQUAL "none") message(FATAL_ERROR "ERROR: Unknown backend type! PARFLOW_ACCELERATOR_BACKEND=${PARFLOW_ACCELERATOR_BACKEND} does not exist!") endif(PARFLOW_ACCELERATOR_BACKEND STREQUAL "none") -# Include Umpire or RMM memory manager for pool allocation +# Include UMPIRE or RMM memory manager for pool allocation if(PARFLOW_HAVE_CUDA OR PARFLOW_HAVE_KOKKOS) - if (DEFINED umpire_ROOT AND DEFINED RMM_ROOT) - message(FATAL_ERROR "ERROR: Cannot have both Umpire and RMM enabled at the same time.") - endif(DEFINED umpire_ROOT AND DEFINED RMM_ROOT) + if (DEFINED UMPIRE_ROOT AND DEFINED RMM_ROOT) + message(FATAL_ERROR "ERROR: Cannot have both UMPIRE and RMM enabled at the same time.") + endif(DEFINED UMPIRE_ROOT AND DEFINED RMM_ROOT) # Determine memory manager and print single status message - if(DEFINED umpire_ROOT) - find_package(umpire REQUIRED) + if(DEFINED UMPIRE_ROOT) + find_package(umpire REQUIRED HINTS ${UMPIRE_ROOT}) set(PARFLOW_HAVE_UMPIRE "yes") - message(STATUS "MEMORY MANAGER: Umpire is sleceted") + message(STATUS "MEMORY MANAGER: UMPIRE is sleceted") - elseif(DEFINED rmm_ROOT) - find_package(rmm REQUIRED) + elseif(DEFINED RMM_ROOT) + find_package(rmm REQUIRED HINTS ${RMM_ROOT}) set(PARFLOW_HAVE_RMM "yes") message(STATUS "MEMORY MANAGER: RMM is sleceted") - else(DEFINED umpire_ROOT) + else(DEFINED UMPIRE_ROOT) message(STATUS "MEMORY MANAGER: NO memoery manager is sleceted") - endif(DEFINED umpire_ROOT) + endif(DEFINED UMPIRE_ROOT) endif(PARFLOW_HAVE_CUDA OR PARFLOW_HAVE_KOKKOS) #----------------------------------------------------------------------------- diff --git a/README-GPU.md b/README-GPU.md index aad504b53..a1f69eb45 100644 --- a/README-GPU.md +++ b/README-GPU.md @@ -9,16 +9,16 @@ Building with CUDA or Kokkos may improve the performance significantly for large ## CMake -Building with GPU acceleration requires a [CUDA](https://docs.nvidia.com/cuda/cuda-installation-guide-linux/index.html) installation, and in case Kokkos is used, [Kokkos](https://github.com/kokkos/kokkos) installation. However, the performance can be further improved by using pool allocation for Unified Memory. Supported memory managers for pool allocation are: [RMM v0.10](https://github.com/rapidsai/rmm/tree/branch-0.10) and [Umpire](https://umpire.readthedocs.io/en/develop/). Note that we are in the process of updating RMM to the newest API, but that should not affect users. Note that you should use only one memory manager, you can't pick both RMM and Umpire. Performance can be improved even more with direct communication between GPUs (requires using a CUDA-Aware MPI library). +Building with GPU acceleration requires a [CUDA](https://docs.nvidia.com/cuda/cuda-installation-guide-linux/index.html) installation, and in case Kokkos is used, [Kokkos](https://github.com/kokkos/kokkos) installation. However, the performance can be further improved by using pool allocation. Supported memory managers for pool allocation are: [RMM](https://github.com/rapidsai/rmm) and [Umpire](https://github.com/llnl/Umpire/tree/develop). Note that you should use only one memory manager, you can't pick both RMM and Umpire. Performance can be improved even more with direct communication between GPUs (requires using a CUDA-Aware MPI library). The GPU acceleration is activated by specifying either *PARFLOW_ACCELERATOR_BACKEND=cuda* option to the CMake, i.e., ```shell cmake ../parflow -DPARFLOW_AMPS_LAYER=mpi1 -DCMAKE_BUILD_TYPE=Release -DPARFLOW_ENABLE_TIMING=TRUE -DPARFLOW_HAVE_CLM=ON -DCMAKE_INSTALL_PREFIX=${PARFLOW_DIR} -DPARFLOW_ACCELERATOR_BACKEND=cuda ``` -or *PARFLOW_ACCELERATOR_BACKEND=kokkos* and *DKokkos_ROOT=/path/to/Kokkos* i.e., +or *PARFLOW_ACCELERATOR_BACKEND=kokkos* and *DKOKKOS_ROOT=/path/to/Kokkos* i.e., ```shell -cmake ../parflow -DPARFLOW_AMPS_LAYER=mpi1 -DCMAKE_BUILD_TYPE=Release -DPARFLOW_ENABLE_TIMING=TRUE -DPARFLOW_HAVE_CLM=ON -DCMAKE_INSTALL_PREFIX=${PARFLOW_DIR} -DPARFLOW_ACCELERATOR_BACKEND=kokkos -DKokkos_ROOT=/path/to/Kokkos +cmake ../parflow -DPARFLOW_AMPS_LAYER=mpi1 -DCMAKE_BUILD_TYPE=Release -DPARFLOW_ENABLE_TIMING=TRUE -DPARFLOW_HAVE_CLM=ON -DCMAKE_INSTALL_PREFIX=${PARFLOW_DIR} -DPARFLOW_ACCELERATOR_BACKEND=kokkos -DKOKKOS_ROOT=/path/to/Kokkos ``` where *DPARFLOW_AMPS_LAYER=mpi1* leverages GPU-based data packing and unpacking. By default, the packed data is copied to a host staging buffer which is then passed for MPI to avoid special requirements for the MPI library. Direct communication between GPUs (with [GPUDirect P2P/RDMA](https://developer.nvidia.com/gpudirect)) can be activated by specifying an environment variable *PARFLOW_USE_GPUDIRECT=1* during runtime in which case the memory copy between CPU and GPU is avoided and a GPU pointer is passed for MPI, but this requires a CUDA-Aware MPI library (support for Unified Memory is not required with the native CUDA backend because the pointers passed to the MPI library point to pinned GPU memory allocations, but is required with the Kokkos backend). @@ -28,17 +28,17 @@ cmake ../parflow -DPARFLOW_AMPS_LAYER=mpi1 -DCMAKE_BUILD_TYPE=Release -DPARFLOW_ ``` or ```shell -cmake ../parflow -DPARFLOW_AMPS_LAYER=mpi1 -DCMAKE_BUILD_TYPE=Release -DPARFLOW_ENABLE_TIMING=TRUE -DPARFLOW_HAVE_CLM=ON -DCMAKE_INSTALL_PREFIX=${PARFLOW_DIR} -DPARFLOW_ACCELERATOR_BACKEND=kokkos -DKokkos_ROOT=/path/to/Kokkos -DRMM_ROOT=/path/to/RMM +cmake ../parflow -DPARFLOW_AMPS_LAYER=mpi1 -DCMAKE_BUILD_TYPE=Release -DPARFLOW_ENABLE_TIMING=TRUE -DPARFLOW_HAVE_CLM=ON -DCMAKE_INSTALL_PREFIX=${PARFLOW_DIR} -DPARFLOW_ACCELERATOR_BACKEND=kokkos -DKOKKOS_ROOT=/path/to/Kokkos -DRMM_ROOT=/path/to/RMM ``` -Similarly, the Umpire library can be activated by specifying the Umpire root directory with -Dumpire_ROOT=/path/to/umpire/root +Similarly, the Umpire library can be activated by specifying the Umpire root directory with -DUMPIRE_ROOT=/path/to/umpire/root ```shell -cmake ../parflow -DPARFLOW_AMPS_LAYER=mpi1 -DCMAKE_BUILD_TYPE=Release -DPARFLOW_ENABLE_TIMING=TRUE -DPARFLOW_HAVE_CLM=ON -DCMAKE_INSTALL_PREFIX=${PARFLOW_DIR} -DPARFLOW_ACCELERATOR_BACKEND=cuda -Dumpire_ROOT=/path/to/umpire +cmake ../parflow -DPARFLOW_AMPS_LAYER=mpi1 -DCMAKE_BUILD_TYPE=Release -DPARFLOW_ENABLE_TIMING=TRUE -DPARFLOW_HAVE_CLM=ON -DCMAKE_INSTALL_PREFIX=${PARFLOW_DIR} -DPARFLOW_ACCELERATOR_BACKEND=cuda -DUMPIRE_ROOT=/path/to/umpire ``` or ```shell -cmake ../parflow -DPARFLOW_AMPS_LAYER=mpi1 -DCMAKE_BUILD_TYPE=Release -DPARFLOW_ENABLE_TIMING=TRUE -DPARFLOW_HAVE_CLM=ON -DCMAKE_INSTALL_PREFIX=${PARFLOW_DIR} -DPARFLOW_ACCELERATOR_BACKEND=kokkos -DKokkos_ROOT=/path/to/Kokkos -Dumpire_ROOT=/path/to/umpire +cmake ../parflow -DPARFLOW_AMPS_LAYER=mpi1 -DCMAKE_BUILD_TYPE=Release -DPARFLOW_ENABLE_TIMING=TRUE -DPARFLOW_HAVE_CLM=ON -DCMAKE_INSTALL_PREFIX=${PARFLOW_DIR} -DPARFLOW_ACCELERATOR_BACKEND=kokkos -DKOKKOS_ROOT=/path/to/Kokkos -DUMPIRE_ROOT=/path/to/umpire ``` Note that on some systems, nvcc cannot locate the MPI include files by default, if this is the case, defining the environment variable CUDAHOSTCXX=mpicxx might help. @@ -47,7 +47,7 @@ Finally, you must make sure you are building the code for the correct GPU archit ```shell -cmake ../parflow -DPARFLOW_AMPS_LAYER=mpi1 -DCMAKE_BUILD_TYPE=Release -DPARFLOW_ENABLE_TIMING=TRUE -DPARFLOW_HAVE_CLM=ON -DCMAKE_INSTALL_PREFIX=${PARFLOW_DIR} -DPARFLOW_ACCELERATOR_BACKEND=cuda -Dumpire_ROOT=/path/to/umpire -DCMAKE_CUDA_ARCHITECTURES=90 +cmake ../parflow -DPARFLOW_AMPS_LAYER=mpi1 -DCMAKE_BUILD_TYPE=Release -DPARFLOW_ENABLE_TIMING=TRUE -DPARFLOW_HAVE_CLM=ON -DCMAKE_INSTALL_PREFIX=${PARFLOW_DIR} -DPARFLOW_ACCELERATOR_BACKEND=cuda -DUMPIRE_ROOT=/path/to/umpire -DCMAKE_CUDA_ARCHITECTURES=90 ``` ## Running Parflow with GPU acceleration diff --git a/docs/user_manual/keys.rst b/docs/user_manual/keys.rst index 4d1f30557..8f99511ca 100644 --- a/docs/user_manual/keys.rst +++ b/docs/user_manual/keys.rst @@ -6337,6 +6337,356 @@ to be active. .Solver.CLM.UseSlopeAspect = True ## Python syntax +.. _CLM Snow Parameterization: + +CLM Snow Parameterization +^^^^^^^^^^^^^^^^^^^^^^^^^ + +These keys control snow physics options for improved snow modeling. +All keys default to backward-compatible behavior when not set. +Note that ``CLM`` must be compiled and linked at runtime for these +options to be active. + +**Rain-Snow Partitioning** + +Standard ``CLM`` uses air temperature to partition precipitation into rain +and snow. However, falling hydrometeors cool via evaporation, making +wet-bulb temperature a better predictor, especially in dry mountain +climates where snow can persist at air temperatures above 0°C. +Reference: :cite:p:`Wang2019`. + +*string* **Solver.CLM.SnowPartition** CLM Selects the method for +partitioning precipitation into rain and snow. The valid types for +this key are **CLM**, **WetbulbThreshold**, **WetbulbLinear**, **Dai**, **Jennings**. + +**CLM**: + Standard air temperature threshold with linear transition (default). + Uses configurable SnowTLow and SnowTHigh thresholds. + +**WetbulbThreshold**: + Sharp threshold at wet-bulb temperature. Better for dry mountain climates. + +**WetbulbLinear**: + Linear transition around wet-bulb temperature threshold. + Uses configurable SnowTransitionWidth for transition zone. + +**Dai**: + Sigmoidal function of air temperature from :cite:p:`Dai2008`. + Uses configurable DaiCoeffA/B/C/D coefficients. + +**Jennings**: + Bivariate logistic regression with air temperature and relative humidity + from :cite:p:`Jennings2018`. Uses configurable JenningsCoeffA/B/G coefficients. + +.. container:: list + + :: + + pfset Solver.CLM.SnowPartition "Dai" ## TCL syntax + .Solver.CLM.SnowPartition = "Dai" ## Python syntax + +*double* **Solver.CLM.SnowTCrit** 2.5 Initial classification threshold +above freezing (K) for determining precipitation type in drv_getforce. +If air temperature exceeds tfrz + SnowTCrit, precipitation is initially +classified as rain. Default 2.5 K matches the hardcoded CLM value. + +.. container:: list + + :: + + pfset Solver.CLM.SnowTCrit 2.5 ## TCL syntax + .Solver.CLM.SnowTCrit = 2.5 ## Python syntax + +*double* **Solver.CLM.SnowTLow** 273.16 CLM method lower temperature +threshold (K) below which all precipitation is snow. Default 273.16 K (freezing). + +.. container:: list + + :: + + pfset Solver.CLM.SnowTLow 273.16 ## TCL syntax + .Solver.CLM.SnowTLow = 273.16 ## Python syntax + +*double* **Solver.CLM.SnowTHigh** 275.16 CLM method upper temperature +threshold (K) above which the liquid fraction reaches maximum (40%). +Default 275.16 K (tfrz + 2). + +.. container:: list + + :: + + pfset Solver.CLM.SnowTHigh 275.16 ## TCL syntax + .Solver.CLM.SnowTHigh = 275.16 ## Python syntax + +*double* **Solver.CLM.SnowTransitionWidth** 1.0 WetbulbLinear method +half-width (K) of transition zone. Transition spans threshold +/- this value. +Default 1.0 K for 2K total range. + +.. container:: list + + :: + + pfset Solver.CLM.SnowTransitionWidth 1.0 ## TCL syntax + .Solver.CLM.SnowTransitionWidth = 1.0 ## Python syntax + +*double* **Solver.CLM.WetbulbThreshold** 274.15 Threshold temperature in +Kelvin for wetbulb partitioning methods. Default 274.15 K (1°C). Only +used when ``Solver.CLM.SnowPartition`` is ``WetbulbThreshold`` or +``WetbulbLinear``. + +.. container:: list + + :: + + pfset Solver.CLM.WetbulbThreshold 274.15 ## TCL syntax + .Solver.CLM.WetbulbThreshold = 274.15 ## Python syntax + +**Dai Coefficients** + +The Dai (2008) method uses a hyperbolic tangent function fitted to 30 years +of global weather station observations. Reference: :cite:p:`Dai2008`. + +Formula: F(%) = a * [tanh(b*(T-c)) - d], where T is in Celsius. +Defaults are from Table 1a (Land, annual) of Dai (2008). + +*double* **Solver.CLM.DaiCoeffA** -48.2292 Dai coefficient a (scaling factor, negative). + +*double* **Solver.CLM.DaiCoeffB** 0.7205 Dai coefficient b (slope at half-frequency T). + +*double* **Solver.CLM.DaiCoeffC** 1.1662 Dai coefficient c (half-frequency temperature in C, +where snow probability is ~50%). + +*double* **Solver.CLM.DaiCoeffD** 1.0223 Dai coefficient d (asymmetry parameter, ~1.0). + +.. container:: list + + :: + + pfset Solver.CLM.DaiCoeffA -48.2292 ## TCL syntax + .Solver.CLM.DaiCoeffA = -48.2292 ## Python syntax + +**Jennings Coefficients** + +The Jennings (2018) method uses bivariate logistic regression: +psnow = 1 / (1 + exp(a + b*T + g*RH)), where T is in Celsius and RH in percent. + +*double* **Solver.CLM.JenningsCoeffA** -10.04 Jennings intercept coefficient. + +*double* **Solver.CLM.JenningsCoeffB** 1.41 Jennings temperature coefficient. + +*double* **Solver.CLM.JenningsCoeffG** 0.09 Jennings relative humidity coefficient. + +.. container:: list + + :: + + pfset Solver.CLM.JenningsCoeffA -10.04 ## TCL syntax + .Solver.CLM.JenningsCoeffA = -10.04 ## Python syntax + +**Thin Snow Damping** + +Thin early-season snowpacks can experience spurious melt due to warm +ground heat flux. This option reduces melt energy for shallow snowpacks +to prevent premature ablation. + +*double* **Solver.CLM.ThinSnowDamping** 1.0 Fraction of melt energy +retained for thin snowpacks. A value of 1.0 means no damping (default, +backward compatible). A value of 0.1 means only 10% of melt energy is +applied (90% reduction). + +.. container:: list + + :: + + pfset Solver.CLM.ThinSnowDamping 0.3 ## TCL syntax + .Solver.CLM.ThinSnowDamping = 0.3 ## Python syntax + +*double* **Solver.CLM.ThinSnowThreshold** 50.0 Snow water equivalent +threshold in mm below which thin snow damping applies. + +.. container:: list + + :: + + pfset Solver.CLM.ThinSnowThreshold 50.0 ## TCL syntax + .Solver.CLM.ThinSnowThreshold = 50.0 ## Python syntax + +**SZA-Based Snow Damping** + +CLM's narrowband snow optical parameters are computed assuming a solar +zenith angle (SZA) of 60 degrees. At higher zenith angles, actual snow +albedo is higher than CLM assumes, meaning less energy should be available +for melt. This is particularly relevant for high-latitude sites during +accumulation season. Reference: :cite:p:`Dang2019`. + +*double* **Solver.CLM.SZASnowDamping** 1.0 Fraction of melt energy +retained at high solar zenith angles. A value of 1.0 means no damping +(default, disabled). A value of 0.8 means 20% energy reduction at high SZA. +Damping varies linearly with cosine of zenith angle between the reference +and minimum thresholds. + +.. container:: list + + :: + + pfset Solver.CLM.SZASnowDamping 0.8 ## TCL syntax + .Solver.CLM.SZASnowDamping = 0.8 ## Python syntax + +*double* **Solver.CLM.SZADampingCoszenRef** 0.5 Reference cosine of solar +zenith angle below which SZA damping applies. Default 0.5 corresponds to +SZA of 60 degrees (matching CLM's assumption for optical parameters). + +.. container:: list + + :: + + pfset Solver.CLM.SZADampingCoszenRef 0.5 ## TCL syntax + .Solver.CLM.SZADampingCoszenRef = 0.5 ## Python syntax + +*double* **Solver.CLM.SZADampingCoszenMin** 0.1 Cosine of solar zenith +angle at which maximum SZA damping applies. Default 0.1 corresponds to +SZA of approximately 84 degrees. Must be less than SZADampingCoszenRef. + +.. container:: list + + :: + + pfset Solver.CLM.SZADampingCoszenMin 0.1 ## TCL syntax + .Solver.CLM.SZADampingCoszenMin = 0.1 ## Python syntax + +Note: Both thin snow damping and SZA damping can be enabled simultaneously. +When both are active, they combine multiplicatively. + +**Albedo Schemes** + +Snow albedo controls the radiation balance and strongly influences melt +timing. Three schemes are available: + +*string* **Solver.CLM.AlbedoScheme** CLM Snow albedo calculation method. +The valid types for this key are **CLM**, **VIC**, **Tarboton**. + +**CLM**: + Age-based exponential decay using the snowage variable (default). + +**VIC**: + Separate decay rates for cold (accumulating) and warm (melting) + conditions based on ground temperature. Reference: :cite:p:`Andreadis2009`. + +**Tarboton**: + Arrhenius temperature-dependent aging where decay accelerates near + the melting point. Reference: :cite:p:`Tarboton1996`. + +.. container:: list + + :: + + pfset Solver.CLM.AlbedoScheme "Tarboton" ## TCL syntax + .Solver.CLM.AlbedoScheme = "Tarboton" ## Python syntax + +*double* **Solver.CLM.AlbedoVisNew** 0.95 Fresh snow visible-band albedo. +Physically ranges 0.85-0.98. + +.. container:: list + + :: + + pfset Solver.CLM.AlbedoVisNew 0.95 ## TCL syntax + .Solver.CLM.AlbedoVisNew = 0.95 ## Python syntax + +*double* **Solver.CLM.AlbedoNirNew** 0.65 Fresh snow near-infrared albedo. +Physically ranges 0.5-0.7. + +.. container:: list + + :: + + pfset Solver.CLM.AlbedoNirNew 0.65 ## TCL syntax + .Solver.CLM.AlbedoNirNew = 0.65 ## Python syntax + +*double* **Solver.CLM.AlbedoMin** 0.4 Minimum snow albedo floor for aged +or dirty snow. + +.. container:: list + + :: + + pfset Solver.CLM.AlbedoMin 0.4 ## TCL syntax + .Solver.CLM.AlbedoMin = 0.4 ## Python syntax + +*double* **Solver.CLM.AlbedoDecayVis** 0.5 Visible albedo decay coefficient +for ``CLM`` and ``Tarboton`` schemes. + +.. container:: list + + :: + + pfset Solver.CLM.AlbedoDecayVis 0.5 ## TCL syntax + .Solver.CLM.AlbedoDecayVis = 0.5 ## Python syntax + +*double* **Solver.CLM.AlbedoDecayNir** 0.2 NIR albedo decay coefficient +for ``CLM`` and ``Tarboton`` schemes. NIR typically decays faster than visible. + +.. container:: list + + :: + + pfset Solver.CLM.AlbedoDecayNir 0.2 ## TCL syntax + .Solver.CLM.AlbedoDecayNir = 0.2 ## Python syntax + +*double* **Solver.CLM.AlbedoAccumA** 0.94 VIC scheme cold-phase +(accumulation) decay base per hour. Should be greater than ``AlbedoThawA``. + +.. container:: list + + :: + + pfset Solver.CLM.AlbedoAccumA 0.94 ## TCL syntax + .Solver.CLM.AlbedoAccumA = 0.94 ## Python syntax + +*double* **Solver.CLM.AlbedoThawA** 0.82 VIC scheme melt-phase decay base +per hour. Should be less than ``AlbedoAccumA`` since melt conditions age +snow faster. + +.. container:: list + + :: + + pfset Solver.CLM.AlbedoThawA 0.82 ## TCL syntax + .Solver.CLM.AlbedoThawA = 0.82 ## Python syntax + +**Fractional Snow Cover** + +The fractional snow covered area (frac_sno) affects surface energy balance +by weighting snow and bare ground contributions. The CLM formulation uses +a tanh-like relationship with snow depth and surface roughness. + +*string* **Solver.CLM.FracSnoScheme** CLM Selects the fractional snow cover +calculation method. Currently only CLM is available; extensible for future +formulations. + +**CLM**: + Standard formulation: frac_sno = snowdp / (10*roughness + snowdp) + +.. container:: list + + :: + + pfset Solver.CLM.FracSnoScheme "CLM" ## TCL syntax + .Solver.CLM.FracSnoScheme = "CLM" ## Python syntax + +*double* **Solver.CLM.FracSnoRoughness** 0.01 Roughness length scale for +fractional snow cover calculation [m]. Default 0.01 m matches CLM's zlnd +parameter for backward compatibility. Larger values reduce snow cover +fraction for a given snow depth. + +.. container:: list + + :: + + pfset Solver.CLM.FracSnoRoughness 0.01 ## TCL syntax + .Solver.CLM.FracSnoRoughness = 0.01 ## Python syntax + + .. _ParFlow NetCDF4 Parallel I/O: ParFlow NetCDF4 Parallel I/O diff --git a/docs/user_manual/refs.bib b/docs/user_manual/refs.bib index aa9377929..e34e9bf7f 100644 --- a/docs/user_manual/refs.bib +++ b/docs/user_manual/refs.bib @@ -256,6 +256,16 @@ @article{Cui14 doi = {10.1016/j.cageo.2014.05.005}, } +@article{Dai2008, +author = {Dai, A.}, +journal = {Geophysical Research Letters}, +pages = {L12802}, +title = {Temperature and pressure dependence of the rain-snow phase transition over land and ocean}, +volume = {35}, +year = {2008}, +doi = {10.1029/2008GL033295}, +} + @article{Dai03, author = {Dai, Y. and X. Zeng and R. E. Dickinson and I. Baker and G. B. Bonan and M. G. Bosilovich and A. S. Denning and P. A. Dirmeyer and P. R., G. Niu and K. W. Oleson and C. A. Schlosser and Z. L. Yang}, journal = {The Bulletin of the American Meteorological Society}, @@ -277,6 +287,16 @@ @article{Danesh-Yazdi2018 url = {https://doi.org/10.1002/hyp.11481}, } +@article{Dang2019, +author = {Dang, C. and Zender, C. S. and Flanner, M. G.}, +journal = {The Cryosphere}, +pages = {2325--2343}, +title = {Intercomparison and improvement of two-stream shortwave radiative transfer schemes in Earth system models for a unified treatment of cryospheric surfaces}, +volume = {13}, +year = {2019}, +doi = {10.5194/tc-13-2325-2019}, +} + @article{DMC10, author = {Daniels, M. H. and Maxwell, R. M. and Chow, F. K.}, journal = {Journal of Hydrologic Engineering}, @@ -557,6 +577,16 @@ @article{Jefferson2017 url = {https://doi.org/10.1175/jhm-d-16-0053.1}, } +@article{Jennings2018, +author = {Jennings, K. S. and Winchell, T. S. and Livneh, B. and Molotch, N. P.}, +journal = {Nature Communications}, +pages = {2831}, +title = {Spatial variation of the rain–snow temperature threshold across the Northern Hemisphere}, +volume = {9}, +year = {2018}, +doi = {10.1038/s41467-018-03629-7}, +} + @article{Jones-Woodward01, author = {Jones, J. E. and Woodward, C. S.}, journal = {Advances in Water Resources}, @@ -1378,3 +1408,36 @@ @misc{endian title = {Endianness --- Wikipedia, The Free Encyclopedia}, url = {http://en.wikipedia.org/wiki/Endianness}, } + +@article{Wang2019, +author = {Wang, Y. and Broxton, P. and Fang, Y. and Behrangi, A. and Barlage, M. and Zeng, X. and Niu, G.}, +journal = {Geophysical Research Letters}, +number = {22}, +pages = {13104--13113}, +title = {A Wet-Bulb Temperature-Based Rain-Snow Partitioning Scheme Improves Snowpack Prediction Over the Drier Western United States}, +volume = {46}, +year = {2019}, +doi = {10.1029/2019GL085722}, +} + +@article{Andreadis2009, +author = {Andreadis, K. M. and Storck, P. and Lettenmaier, D. P.}, +journal = {Journal of Hydrometeorology}, +number = {6}, +pages = {1379--1395}, +title = {Modeling Snow Accumulation and Ablation Processes in Forested Environments}, +volume = {10}, +year = {2009}, +doi = {10.1175/2009JHM1083.1}, +} + +@article{Tarboton1996, +author = {Tarboton, D. G. and Luce, C. H.}, +journal = {Water Resources Research}, +number = {4}, +pages = {1195--1207}, +title = {Utah Energy Balance Snow Accumulation and Melt Model (UEB)}, +volume = {32}, +year = {1996}, +doi = {10.1029/96WR01049}, +} diff --git a/pf-keys/definitions/solver.yaml b/pf-keys/definitions/solver.yaml index 25a029370..d8cf67bf3 100644 --- a/pf-keys/definitions/solver.yaml +++ b/pf-keys/definitions/solver.yaml @@ -494,6 +494,340 @@ Solver: BoolDomain: RequiresModule: CLM + # ----------------------------------------------------------------------------- + # CLM Snow Parameterization Keys + # ----------------------------------------------------------------------------- + + SnowPartition: + help: > + [Type: string] Selects the rain-snow partitioning method. CLM uses air + temperature threshold, while wetbulb methods account for evaporative + cooling of falling hydrometeors. Dai (2008) uses a sigmoidal function + of air temperature. Jennings (2018) uses bivariate logistic regression + with air temperature and relative humidity. + References: Wang et al. (2019) GRL, Dai (2008) JAMC, Jennings et al. (2018) Nat Commun. + default: CLM + domains: + EnumDomain: + enum_list: + - CLM + - WetbulbThreshold + - WetbulbLinear + - Dai + - Jennings + RequiresModule: CLM + + SnowTCrit: + help: > + [Type: double] Initial classification threshold above tfrz (K) for determining + precipitation type in drv_getforce. Default 2.5 K matches hardcoded CLM value. + If T > tfrz + SnowTCrit, precipitation is initially classified as rain. + Further partitioning occurs in clm_hydro_canopy based on SnowPartition method. + default: 2.5 + domains: + DoubleValue: + min_value: 0.0 + max_value: 5.0 + RequiresModule: CLM + + SnowTLow: + help: > + [Type: double] CLM method lower temperature threshold (K) below which all + precipitation is snow. Default 273.16 K (freezing point). + default: 273.16 + domains: + DoubleValue: + min_value: 268.16 + max_value: 278.16 + RequiresModule: CLM + + SnowTHigh: + help: > + [Type: double] CLM method upper temperature threshold (K) above which + liquid fraction reaches maximum (40%). Default 275.16 K (tfrz + 2). + default: 275.16 + domains: + DoubleValue: + min_value: 270.16 + max_value: 280.16 + RequiresModule: CLM + + SnowTransitionWidth: + help: > + [Type: double] WetbulbLinear method half-width (K) of transition zone. + Transition spans threshold +/- this value. Default 1.0 K for 2K total range. + default: 1.0 + domains: + DoubleValue: + min_value: 0.1 + max_value: 5.0 + RequiresModule: CLM + + DaiCoeffA: + help: > + [Type: double] Dai (2008) coefficient a (scaling factor, negative). + F(%) = a * [tanh(b*(T-c)) - d]. Default -48.2292 from Table 1a (Land, ANN). + Reference: Dai (2008) GRL doi:10.1029/2008GL033295 + default: -48.2292 + domains: + DoubleValue: + min_value: -100.0 + max_value: 0.0 + RequiresModule: CLM + + DaiCoeffB: + help: > + [Type: double] Dai (2008) coefficient b (slope at half-frequency T). + F(%) = a * [tanh(b*(T-c)) - d]. Default 0.7205 from Table 1a (Land, ANN). + default: 0.7205 + domains: + DoubleValue: + min_value: 0.1 + max_value: 2.0 + RequiresModule: CLM + + DaiCoeffC: + help: > + [Type: double] Dai (2008) coefficient c (half-frequency temperature in C). + F(%) = a * [tanh(b*(T-c)) - d]. Default 1.1662 C from Table 1a (Land, ANN). + This is the temperature at which snow probability is approximately 50%. + default: 1.1662 + domains: + DoubleValue: + min_value: -2.0 + max_value: 5.0 + RequiresModule: CLM + + DaiCoeffD: + help: > + [Type: double] Dai (2008) coefficient d (asymmetry parameter). + F(%) = a * [tanh(b*(T-c)) - d]. Default 1.0223 from Table 1a (Land, ANN). + Values close to 1.0 indicate symmetric transition. + default: 1.0223 + domains: + DoubleValue: + min_value: 0.5 + max_value: 1.5 + RequiresModule: CLM + + JenningsCoeffA: + help: > + [Type: double] Jennings et al. (2018) bivariate logistic intercept. + psnow = 1 / (1 + exp(a + b*T + g*RH)). Default from Jennings Table 1. + default: -10.04 + domains: + DoubleValue: + min_value: -20.0 + max_value: 0.0 + RequiresModule: CLM + + JenningsCoeffB: + help: > + [Type: double] Jennings et al. (2018) bivariate logistic temperature coefficient. + psnow = 1 / (1 + exp(a + b*T + g*RH)). T in Celsius. Default from Jennings Table 1. + default: 1.41 + domains: + DoubleValue: + min_value: 0.1 + max_value: 3.0 + RequiresModule: CLM + + JenningsCoeffG: + help: > + [Type: double] Jennings et al. (2018) bivariate logistic relative humidity coefficient. + psnow = 1 / (1 + exp(a + b*T + g*RH)). RH in percent. Default from Jennings Table 1. + default: 0.09 + domains: + DoubleValue: + min_value: 0.01 + max_value: 0.5 + RequiresModule: CLM + + WetbulbThreshold: + help: > + [Type: double] Threshold temperature (K) for wetbulb rain-snow partitioning. + Default 274.15 K (1 degree C). Only used when SnowPartition is WetbulbThreshold + or WetbulbLinear. + default: 274.15 + domains: + DoubleValue: + min_value: 268.15 + max_value: 278.15 + RequiresModule: CLM + + ThinSnowDamping: + help: > + [Type: double] Fraction of melt energy retained for thin snowpacks. + 1.0 = no damping (default, backward compatible), 0.1 = 90% energy reduction. + Prevents spurious early-season melt from warm ground heat flux. + default: 1.0 + domains: + DoubleValue: + min_value: 0.0 + max_value: 1.0 + RequiresModule: CLM + + ThinSnowThreshold: + help: > + [Type: double] SWE threshold (mm) below which thin snow damping applies. + default: 50.0 + domains: + DoubleValue: + min_value: 1.0 + max_value: 500.0 + RequiresModule: CLM + + SZASnowDamping: + help: > + [Type: double] Fraction of melt energy retained at high solar zenith angles. + 1.0 = no damping (default, disabled), 0.8 = 20% energy reduction at high SZA. + Compensates for CLM underestimating albedo at SZA > 60 degrees. + Physical basis: CLM narrowband parameters assume SZA=60 deg. + Reference: Dang et al. (2019) The Cryosphere. + default: 1.0 + domains: + DoubleValue: + min_value: 0.0 + max_value: 1.0 + RequiresModule: CLM + + SZADampingCoszenRef: + help: > + [Type: double] Reference cosine of solar zenith angle below which SZA damping + applies. Default 0.5 corresponds to SZA=60 degrees (CLM assumption). + Damping scales linearly from 1.0 at this coszen to SZASnowDamping at CoszenMin. + default: 0.5 + domains: + DoubleValue: + min_value: 0.2 + max_value: 1.0 + RequiresModule: CLM + + SZADampingCoszenMin: + help: > + [Type: double] Cosine of solar zenith angle at which maximum SZA damping applies. + Default 0.1 corresponds to SZA of approximately 84 degrees. + Must be less than SZADampingCoszenRef. + default: 0.1 + domains: + DoubleValue: + min_value: 0.0 + max_value: 0.5 + RequiresModule: CLM + + AlbedoScheme: + help: > + [Type: string] Snow albedo calculation method. CLM uses age-based decay, + VIC uses separate cold/melt decay rates based on ground temperature, + Tarboton uses Arrhenius temperature-dependent aging. + References: Andreadis et al. (2009), Tarboton and Luce (1996). + default: CLM + domains: + EnumDomain: + enum_list: + - CLM + - VIC + - Tarboton + RequiresModule: CLM + + AlbedoVisNew: + help: > + [Type: double] Fresh snow visible-band albedo. Physically ranges 0.85-0.98. + default: 0.95 + domains: + DoubleValue: + min_value: 0.80 + max_value: 0.99 + RequiresModule: CLM + + AlbedoNirNew: + help: > + [Type: double] Fresh snow near-infrared albedo. Physically ranges 0.5-0.7. + default: 0.65 + domains: + DoubleValue: + min_value: 0.40 + max_value: 0.85 + RequiresModule: CLM + + AlbedoMin: + help: > + [Type: double] Minimum snow albedo floor for aged or dirty snow. + default: 0.4 + domains: + DoubleValue: + min_value: 0.2 + max_value: 0.7 + RequiresModule: CLM + + AlbedoDecayVis: + help: > + [Type: double] Visible albedo decay coefficient for CLM and Tarboton schemes. + default: 0.5 + domains: + DoubleValue: + min_value: 0.1 + max_value: 2.0 + RequiresModule: CLM + + AlbedoDecayNir: + help: > + [Type: double] NIR albedo decay coefficient for CLM and Tarboton schemes. + NIR typically decays faster than visible. + default: 0.2 + domains: + DoubleValue: + min_value: 0.05 + max_value: 1.0 + RequiresModule: CLM + + AlbedoAccumA: + help: > + [Type: double] VIC scheme cold-phase (accumulation) decay base per hour. + Should be greater than AlbedoThawA. + default: 0.94 + domains: + DoubleValue: + min_value: 0.80 + max_value: 0.99 + RequiresModule: CLM + + AlbedoThawA: + help: > + [Type: double] VIC scheme melt-phase decay base per hour. + Should be less than AlbedoAccumA since melt conditions age snow faster. + default: 0.82 + domains: + DoubleValue: + min_value: 0.70 + max_value: 0.95 + RequiresModule: CLM + + FracSnoScheme: + help: > + [Type: string] Selects the fractional snow covered area (frac_sno) calculation method. + CLM uses the standard formulation: frac_sno = snowdp / (10*roughness + snowdp). + Additional formulations may be added in future versions. + default: CLM + domains: + EnumDomain: + enum_list: + - CLM + RequiresModule: CLM + + FracSnoRoughness: + help: > + [Type: double] Roughness length scale for fractional snow cover calculation [m]. + Default 0.01 m matches CLM's zlnd parameter for backward compatibility. + Larger values reduce snow cover fraction for a given snow depth. + Formula: frac_sno = snowdp / (10*FracSnoRoughness + snowdp) + default: 0.01 + domains: + DoubleValue: + min_value: 0.001 + max_value: 0.1 + RequiresModule: CLM + # ----------------------------------------------------------------------------- # CLM input variables to write drv_clmin.dat (Input) and drv_vegp.dat (VegParams) files # ----------------------------------------------------------------------------- diff --git a/pfsimulator/clm/clm.F90 b/pfsimulator/clm/clm.F90 index b2c2cb994..4ae52c07c 100644 --- a/pfsimulator/clm/clm.F90 +++ b/pfsimulator/clm/clm.F90 @@ -10,7 +10,14 @@ subroutine clm_lsm(pressure,saturation,evap_trans,top,bottom,porosity,pf_dz_mult write_CLM_binary,slope_accounting_CLM,beta_typepf,veg_water_stress_typepf,wilting_pointpf,field_capacitypf, & res_satpf,irr_typepf, irr_cyclepf, irr_ratepf, irr_startpf, irr_stoppf, irr_thresholdpf, & qirr_pf,qirr_inst_pf,irr_flag_pf,irr_thresholdtypepf,soi_z,clm_next,clm_write_logs, & -clm_last_rst,clm_daily_rst,rz_water_stress_typepf, pf_nlevsoi, pf_nlevlak) +clm_last_rst,clm_daily_rst,rz_water_stress_typepf, pf_nlevsoi, pf_nlevlak, & +snow_partition_typepf,tw_thresholdpf,thin_snow_dampingpf,thin_snow_thresholdpf, & +snow_tcritpf,snow_t_lowpf,snow_t_highpf,snow_transition_widthpf, & +dai_apf,dai_bpf,dai_cpf,dai_dpf,jennings_apf,jennings_bpf,jennings_gpf, & +sza_snow_dampingpf,sza_damping_coszen_refpf,sza_damping_coszen_minpf, & +albedo_schemepf,albedo_vis_newpf,albedo_nir_newpf,albedo_minpf, & +albedo_decay_vispf,albedo_decay_nirpf,albedo_accum_apf,albedo_thaw_apf, & +frac_sno_typepf,frac_sno_roughnesspf) !========================================================================= ! @@ -147,6 +154,42 @@ subroutine clm_lsm(pressure,saturation,evap_trans,top,bottom,porosity,pf_dz_mult real(r8) :: irr_thresholdpf ! irrigation threshold criteria for deficit cycle (units of soil moisture content) integer :: irr_thresholdtypepf ! irrigation threshold criteria type -- top layer, bottom layer, column avg + ! snow parameterization keys @RMM 2025 + integer :: snow_partition_typepf ! rain-snow partition: 0=CLM, 1=wb thresh, 2=wb lin, 3=Dai, 4=Jennings + real(r8) :: tw_thresholdpf ! wetbulb temperature threshold for snow [K] + real(r8) :: thin_snow_dampingpf ! thin snow energy damping factor [0-1] + real(r8) :: thin_snow_thresholdpf ! SWE threshold for damping [kg/m2] + real(r8) :: snow_tcritpf ! initial T classification threshold above tfrz [K], default 2.5 + real(r8) :: snow_t_lowpf ! CLM method lower T threshold [K], default 273.16 + real(r8) :: snow_t_highpf ! CLM method upper T threshold [K], default 275.16 + real(r8) :: snow_transition_widthpf ! WetbulbLinear half-width [K], default 1.0 + real(r8) :: dai_apf ! Dai (2008) coefficient a, default -48.2292 + real(r8) :: dai_bpf ! Dai (2008) coefficient b, default 0.7205 + real(r8) :: dai_cpf ! Dai (2008) coefficient c, default 1.1662 + real(r8) :: dai_dpf ! Dai (2008) coefficient d, default 1.0223 + real(r8) :: jennings_apf ! Jennings (2018) intercept, default -10.04 + real(r8) :: jennings_bpf ! Jennings (2018) T coefficient, default 1.41 + real(r8) :: jennings_gpf ! Jennings (2018) RH coefficient, default 0.09 + + ! SZA snow damping keys @RMM 2025 + real(r8) :: sza_snow_dampingpf ! SZA damping factor [0-1], 1.0=disabled + real(r8) :: sza_damping_coszen_refpf ! reference coszen for damping onset + real(r8) :: sza_damping_coszen_minpf ! coszen at max damping + + ! snow albedo parameterization keys @RMM 2025 + integer :: albedo_schemepf ! albedo scheme: 0=CLM, 1=VIC, 2=Tarboton + real(r8) :: albedo_vis_newpf ! fresh snow VIS albedo [0-1] + real(r8) :: albedo_nir_newpf ! fresh snow NIR albedo [0-1] + real(r8) :: albedo_minpf ! minimum albedo floor [0-1] + real(r8) :: albedo_decay_vispf ! VIS decay coefficient [0-1] + real(r8) :: albedo_decay_nirpf ! NIR decay coefficient [0-1] + real(r8) :: albedo_accum_apf ! VIC cold-phase decay base + real(r8) :: albedo_thaw_apf ! VIC melt-phase decay base + + ! frac_sno parameterization keys @RMM 2025 + integer :: frac_sno_typepf ! frac_sno scheme: 0=CLM (default), others TBD + real(r8) :: frac_sno_roughnesspf ! roughness length for frac_sno [m] + ! local indices & counters integer :: i,j,k,k1,j1,l1 ! indices for local looping integer :: bj,bl ! indices for local looping !BH @@ -487,6 +530,43 @@ subroutine clm_lsm(pressure,saturation,evap_trans,top,bottom,porosity,pf_dz_mult clm(t)%field_capacity = field_capacitypf clm(t)%res_sat = res_satpf + ! for snow parameterization @RMM 2025 + clm(t)%snow_partition_type = snow_partition_typepf + clm(t)%tw_threshold = tw_thresholdpf + clm(t)%thin_snow_damping = thin_snow_dampingpf + clm(t)%thin_snow_threshold = thin_snow_thresholdpf + clm(t)%snow_tcrit = snow_tcritpf + clm(t)%snow_t_low = snow_t_lowpf + clm(t)%snow_t_high = snow_t_highpf + clm(t)%snow_transition_width = snow_transition_widthpf + clm(t)%dai_a = dai_apf + clm(t)%dai_b = dai_bpf + clm(t)%dai_c = dai_cpf + clm(t)%dai_d = dai_dpf + clm(t)%jennings_a = jennings_apf + clm(t)%jennings_b = jennings_bpf + clm(t)%jennings_g = jennings_gpf + + ! for SZA snow damping @RMM 2025 + clm(t)%sza_snow_damping = sza_snow_dampingpf + clm(t)%sza_damping_coszen_ref = sza_damping_coszen_refpf + clm(t)%sza_damping_coszen_min = sza_damping_coszen_minpf + clm(t)%coszen = 0.5d0 ! initialize to reference value + + ! for snow albedo parameterization @RMM 2025 + clm(t)%albedo_scheme = albedo_schemepf + clm(t)%albedo_vis_new = albedo_vis_newpf + clm(t)%albedo_nir_new = albedo_nir_newpf + clm(t)%albedo_min = albedo_minpf + clm(t)%albedo_decay_vis = albedo_decay_vispf + clm(t)%albedo_decay_nir = albedo_decay_nirpf + clm(t)%albedo_accum_a = albedo_accum_apf + clm(t)%albedo_thaw_a = albedo_thaw_apf + + ! for frac_sno parameterization @RMM 2025 + clm(t)%frac_sno_type = frac_sno_typepf + clm(t)%frac_sno_roughness = frac_sno_roughnesspf + ! for irrigation clm(t)%irr_type = irr_typepf clm(t)%irr_cycle = irr_cyclepf diff --git a/pfsimulator/clm/clm_dynvegpar.F90 b/pfsimulator/clm/clm_dynvegpar.F90 index 90fce856b..e713b4e09 100644 --- a/pfsimulator/clm/clm_dynvegpar.F90 +++ b/pfsimulator/clm/clm_dynvegpar.F90 @@ -79,7 +79,18 @@ subroutine clm_dynvegpar (clm,clm_forc_veg) endif ! Fraction of soil covered by snow +! @RMM 2025: Added configurable frac_sno schemes and adjustable roughness parameter - clm%frac_sno = clm%snowdp/(10.*clm%zlnd + clm%snowdp) + select case (clm%frac_sno_type) + + case (0) ! CLM default + ! Use configurable roughness parameter (defaults to zlnd for backward compatibility) + clm%frac_sno = clm%snowdp / (10.0d0 * clm%frac_sno_roughness + clm%snowdp) + + case default ! Future formulations TBD + ! Default to CLM formulation + clm%frac_sno = clm%snowdp / (10.0d0 * clm%frac_sno_roughness + clm%snowdp) + + end select end subroutine clm_dynvegpar diff --git a/pfsimulator/clm/clm_hydro_canopy.F90 b/pfsimulator/clm/clm_hydro_canopy.F90 index a065b1cd8..a5b2220ad 100644 --- a/pfsimulator/clm/clm_hydro_canopy.F90 +++ b/pfsimulator/clm/clm_hydro_canopy.F90 @@ -32,6 +32,11 @@ subroutine clm_hydro_canopy (clm) use clm_varcon, only : tfrz, istice, istwet, istsoil implicit none + ! Parameters for saturation vapor pressure calculation + real(r8), parameter :: es_a = 611.2d0 ! [Pa] reference saturation vapor pressure + real(r8), parameter :: es_b = 17.67d0 ! coefficient for Clausius-Clapeyron + real(r8), parameter :: es_c = 243.5d0 ! [C] coefficient for Clausius-Clapeyron + !=== Arguments =========================================================== type (clm1d), intent(inout) :: clm !CLM 1-D Module @@ -56,6 +61,17 @@ subroutine clm_hydro_canopy (clm) real(r8) :: & dewmxi ! inverse of maximum allowed dew [1/mm] + ! Variables for wetbulb rain-snow partitioning @RMM 2025 + real(r8) :: t_c ! air temperature in Celsius + real(r8) :: t_wb ! wetbulb temperature in Celsius + real(r8) :: t_wb_k ! wetbulb temperature in Kelvin + real(r8) :: rh_pct ! relative humidity in percent + real(r8) :: e_sat ! saturation vapor pressure [Pa] + real(r8) :: e_act ! actual vapor pressure [Pa] + real(r8) :: q_sat ! saturation specific humidity [kg/kg] + real(r8) :: psnow ! probability/fraction of snow [-] + real(r8) :: exponent ! exponent for logistic function + !=== End Variable List =================================================== qflx_candrip = 0. @@ -144,22 +160,117 @@ subroutine clm_hydro_canopy (clm) clm%qflx_prec_grnd = qflx_through + qflx_candrip endif - ! 1.3 The percentage of liquid water by mass, which is arbitrarily set to + ! 1.3 The percentage of liquid water by mass, which is arbitrarily set to ! vary linearly with air temp, from 0% at 273.16 to 40% max at 275.16. + ! @RMM 2025: Only use itypprc shortcut for CLM default method (type 0). + ! Wetbulb methods (type 1,2) need to run regardless of itypprc to catch + ! cases where air temp is warm but wetbulb is cold (humid conditions). - if (clm%itypprc <= 1) then + if (clm%itypprc <= 1 .and. clm%snow_partition_type == 0) then flfall = 1. ! fraction of liquid water within falling precip. clm%qflx_snow_grnd = 0. ! ice onto ground (mm/s) clm%qflx_rain_grnd = clm%qflx_prec_grnd ! liquid water onto ground (mm/s) dz_snowf = 0. ! rate of snowfall, snow depth/s (m/s) else - if (clm%forc_t <= tfrz) then - flfall = 0. - else if (clm%forc_t <= tfrz+2.) then - flfall = -54.632 + 0.2*clm%forc_t - else - flfall = 0.4 - endif + ! @RMM 2025: Select rain-snow partitioning method + ! snow_partition_type: 0=CLM linear, 1=wetbulb thresh, 2=wetbulb linear, 3=Dai, 4=Jennings + select case (clm%snow_partition_type) + + case (1, 2) ! Wetbulb-based methods + ! Calculate wetbulb temperature using Stull (2011) psychrometric approximation + ! First convert air temp to Celsius + t_c = clm%forc_t - tfrz + + ! Calculate saturation vapor pressure using Clausius-Clapeyron + e_sat = es_a * exp(es_b * t_c / (t_c + es_c)) + + ! Calculate saturation specific humidity + q_sat = 0.622d0 * e_sat / (clm%forc_pbot - 0.378d0 * e_sat) + + ! Calculate relative humidity from specific humidity + ! Avoid division by zero and cap at 100% + if (q_sat > 0.0d0) then + rh_pct = min(100.0d0, max(0.0d0, 100.0d0 * clm%forc_q / q_sat)) + else + rh_pct = 100.0d0 + endif + + ! Stull (2011) wet-bulb temperature approximation (result in Celsius) + ! Valid for RH >= 5% and T between -20C and 50C + t_wb = t_c * atan(0.151977d0 * sqrt(rh_pct + 8.313659d0)) & + + atan(t_c + rh_pct) & + - atan(rh_pct - 1.676331d0) & + + 0.00391838d0 * (rh_pct**1.5d0) * atan(0.023101d0 * rh_pct) & + - 4.686035d0 + + ! Convert wetbulb to Kelvin for comparison + t_wb_k = t_wb + tfrz + + if (clm%snow_partition_type == 1) then + ! Wetbulb threshold method: all snow below threshold, all rain above + if (t_wb_k <= clm%tw_threshold) then + flfall = 0.0d0 ! all snow + else + flfall = 1.0d0 ! all rain + endif + else ! case 2: wetbulb linear + ! Linear transition over configurable range centered on threshold + if (t_wb_k <= clm%tw_threshold - clm%snow_transition_width) then + flfall = 0.0d0 + else if (t_wb_k >= clm%tw_threshold + clm%snow_transition_width) then + flfall = 1.0d0 + else + flfall = (t_wb_k - (clm%tw_threshold - clm%snow_transition_width)) / & + (2.0d0 * clm%snow_transition_width) + endif + endif + + case (3) ! Dai (2008) sigmoidal method + ! F(%) = a * [tanh(b*(T-c)) - d], converted to fraction by /100 + ! T in Celsius, coefficients from Table 1a (Land, ANN) + ! Reference: Dai (2008) GRL doi:10.1029/2008GL033295 + t_c = clm%forc_t - tfrz + psnow = (clm%dai_a / 100.0d0) * & + (tanh(clm%dai_b * (t_c - clm%dai_c)) - clm%dai_d) + psnow = max(0.0d0, min(1.0d0, psnow)) + flfall = 1.0d0 - psnow + + case (4) ! Jennings et al. (2018) bivariate logistic method + ! psnow = 1 / (1 + exp(a + b*T + g*RH)) + ! T in Celsius, RH in percent + ! Reference: Jennings et al. (2018) Nat Commun doi:10.1038/s41467-018-03629-7 + t_c = clm%forc_t - tfrz + + ! Calculate saturation vapor pressure using Clausius-Clapeyron + e_sat = es_a * exp(es_b * t_c / (t_c + es_c)) + + ! Calculate saturation specific humidity + q_sat = 0.622d0 * e_sat / (clm%forc_pbot - 0.378d0 * e_sat) + + ! Calculate relative humidity from specific humidity + if (q_sat > 0.0d0) then + rh_pct = min(100.0d0, max(0.0d0, 100.0d0 * clm%forc_q / q_sat)) + else + rh_pct = 100.0d0 + endif + + ! Bivariate logistic regression + exponent = clm%jennings_a + clm%jennings_b * t_c + clm%jennings_g * rh_pct + psnow = 1.0d0 / (1.0d0 + exp(exponent)) + flfall = 1.0d0 - psnow + + case default ! Case 0: CLM default linear (air temperature) with configurable thresholds + ! Now uses configurable snow_t_low and snow_t_high instead of hardcoded values + if (clm%forc_t <= clm%snow_t_low) then + flfall = 0.0d0 + else if (clm%forc_t >= clm%snow_t_high) then + flfall = 0.4d0 + else + flfall = 0.4d0 * (clm%forc_t - clm%snow_t_low) / & + (clm%snow_t_high - clm%snow_t_low) + endif + + end select ! Use Alta relationship, Anderson(1976); LaChapelle(1961), ! U.S.Department of Agriculture Forest Service, Project F, diff --git a/pfsimulator/clm/clm_main.F90 b/pfsimulator/clm/clm_main.F90 index 30205855c..cc66ec6c0 100644 --- a/pfsimulator/clm/clm_main.F90 +++ b/pfsimulator/clm/clm_main.F90 @@ -184,6 +184,9 @@ subroutine clm_main (clm,day,gmt,clm_forc_veg) call clm_coszen (clm,day,coszen) + ! @RMM 2025: Store coszen for use in meltfreeze SZA damping + clm%coszen = coszen + call clm_surfalb (clm,coszen) ! ----------------------------------------------------------------- diff --git a/pfsimulator/clm/clm_meltfreeze.F90 b/pfsimulator/clm/clm_meltfreeze.F90 index 691a3c180..9f95728f6 100644 --- a/pfsimulator/clm/clm_meltfreeze.F90 +++ b/pfsimulator/clm/clm_meltfreeze.F90 @@ -62,7 +62,12 @@ subroutine clm_meltfreeze (fact, brr, hs, dhsdT, & temp1 ! temporary variables [kg/m2] real(r8), dimension(clm%snl+1 : nlevsoi) :: wmass0, wice0, wliq0 - real(r8) propor,tinc + real(r8) propor,tinc + + ! @RMM 2025: Combined thin snow + SZA damping variables + real(r8) :: damping_factor ! combined damping factor + real(r8) :: depth_damping_factor ! depth-based (thin snow) damping component + real(r8) :: sza_damping_factor ! SZA-based damping component !=== End Variable List =================================================== @@ -119,6 +124,55 @@ subroutine clm_meltfreeze (fact, brr, hs, dhsdT, & endif enddo +! @RMM 2025: Apply combined thin-snow and SZA-based damping to snow layer melt energy +! Both effects can be enabled independently and combine multiplicatively: +! - Thin snow damping: reduces melt when SWE is below threshold +! - SZA damping: reduces melt at high solar zenith angles where CLM underestimates albedo + do j = clm%snl+1, 0 ! Only snow layers (indices <= 0) + if (clm%imelt(j) == 1 .and. hm(j) > 0.0d0) then ! Only for melting + + ! Initialize both factors to 1.0 (no damping) + depth_damping_factor = 1.0d0 + sza_damping_factor = 1.0d0 + + ! Calculate depth-based (thin snow) damping + ! Linear interpolation from damping at SWE=0 to 1.0 at threshold + if (clm%thin_snow_damping > 0.0d0 .and. clm%thin_snow_damping < 1.0d0) then + if (clm%h2osno <= 0.0d0) then + depth_damping_factor = clm%thin_snow_damping + else if (clm%h2osno >= clm%thin_snow_threshold) then + depth_damping_factor = 1.0d0 + else + depth_damping_factor = clm%thin_snow_damping + & + (1.0d0 - clm%thin_snow_damping) * & + (clm%h2osno / clm%thin_snow_threshold) + endif + endif + + ! Calculate SZA-based damping + ! At coszen >= coszen_ref (low SZA, e.g. 60 deg): no damping + ! At coszen <= coszen_min (high SZA, e.g. 84 deg): maximum damping + ! Linear interpolation between bounds + if (clm%sza_snow_damping > 0.0d0 .and. clm%sza_snow_damping < 1.0d0) then + if (clm%coszen >= clm%sza_damping_coszen_ref) then + sza_damping_factor = 1.0d0 + else if (clm%coszen <= clm%sza_damping_coszen_min) then + sza_damping_factor = clm%sza_snow_damping + else + sza_damping_factor = clm%sza_snow_damping + & + (1.0d0 - clm%sza_snow_damping) * & + (clm%coszen - clm%sza_damping_coszen_min) / & + (clm%sza_damping_coszen_ref - clm%sza_damping_coszen_min) + endif + endif + + ! Combine multiplicatively and apply to melt energy + damping_factor = depth_damping_factor * sza_damping_factor + hm(j) = hm(j) * damping_factor + + endif + enddo + ! These two errors were checked carefully. They result from the the computed error ! of "Tridiagonal-Matrix" in subroutine "thermal". diff --git a/pfsimulator/clm/clm_snowalb.F90 b/pfsimulator/clm/clm_snowalb.F90 index b441a1407..6b7503a4e 100644 --- a/pfsimulator/clm/clm_snowalb.F90 +++ b/pfsimulator/clm/clm_snowalb.F90 @@ -2,20 +2,27 @@ subroutine clm_snowalb (clm, coszen, nband, ind, alb) -!----------------------------------------------------------------------- -! -! Purpose: -! Determine snow albedos -! +!----------------------------------------------------------------------- +! +! Purpose: +! Determine snow albedos with configurable parameterization schemes +! +! Schemes available: +! 0 = CLM (default) - age-based decay with configurable parameters +! 1 = VIC - dual cold/warm decay rates based on ground temperature +! 2 = Tarboton - Arrhenius temperature-dependent aging +! !----------------------------------------------------------------------- ! $Id: clm_snowalb.F90,v 1.1.1.1 2006/02/14 23:05:52 kollet Exp $ +! Modified @RMM 2025 - configurable albedo parameters and alternative schemes !----------------------------------------------------------------------- use precision use clmtype + use clm_varcon, only : tfrz implicit none -! ------------------------ arguments ------------------------------ +! ------------------------ arguments ------------------------------ type (clm1d), intent(inout) :: clm !CLM 1-D Module real(r8) , intent(in) :: coszen !cosine solar zenith angle for next time step integer , intent(in) :: nband !number of solar radiation waveband classes @@ -27,13 +34,16 @@ subroutine clm_snowalb (clm, coszen, nband, ind, alb) integer :: ib !waveband class - real(r8) :: snal0 = 0.95 !vis albedo of new snow for sza<60 - real(r8) :: snal1 = 0.65 !nir albedo of new snow for sza<60 - real(r8) :: conn = 0.5 !constant for visible snow alb calculation [-] - real(r8) :: cons = 0.2 !constant (=0.2) for nir snow albedo calculation [-] real(r8) :: sl = 2.0 !factor that helps control alb zenith dependence [-] + real(r8) :: snal0 !vis albedo of new snow (from clm state) + real(r8) :: snal1 !nir albedo of new snow (from clm state) + real(r8) :: conn !constant for visible snow alb calculation (from clm state) + real(r8) :: cons !constant for nir snow albedo calculation (from clm state) + real(r8) :: age !factor to reduce visible snow alb due to snow age [-] + real(r8) :: age_days !snow age in days for VIC/Tarboton schemes + real(r8) :: decay_factor !temperature-dependent decay factor for Tarboton scheme real(r8) :: albs !temporary vis snow albedo real(r8) :: albl !temporary nir snow albedo real(r8) :: cff !snow alb correction factor for zenith angle > 60 [-] @@ -41,6 +51,12 @@ subroutine clm_snowalb (clm, coszen, nband, ind, alb) ! ----------------------------------------------------------------- +! Use configurable parameters from clm state + snal0 = clm%albedo_vis_new ! default 0.95 + snal1 = clm%albedo_nir_new ! default 0.65 + conn = clm%albedo_decay_vis ! default 0.5 + cons = clm%albedo_decay_nir ! default 0.2 + ! zero albedos do ib = 1, nband @@ -48,27 +64,93 @@ subroutine clm_snowalb (clm, coszen, nband, ind, alb) end do ! ========================================================================= -! CLM Albedo for snow cover. -! Snow albedo depends on snow-age, zenith angle, and thickness of snow age -! gives reduction of visible radiation +! Snow albedo calculation based on selected scheme ! ========================================================================= -! Correction for snow age + select case (clm%albedo_scheme) + + case (0) ! CLM default - age-based decay + ! ===================================================================== + ! CLM Albedo for snow cover. + ! Snow albedo depends on snow-age, zenith angle, and thickness of snow age + ! gives reduction of visible radiation + ! ===================================================================== + + ! Correction for snow age + age = 1._r8 - 1._r8/(1._r8 + clm%snowage) + albs = snal0*(1._r8 - cons*age) + albl = snal1*(1._r8 - conn*age) + + case (1) ! VIC - dual cold/warm decay rates + ! ===================================================================== + ! VIC-style albedo decay with different rates for cold vs warm conditions + ! Cold (accumulating): slower decay, albedo_accum_a base + ! Warm (melting): faster decay, albedo_thaw_a base + ! Reference: Andreadis et al. (2009), VIC snow model documentation + ! ===================================================================== + + ! Convert snow age to days (snowage is in seconds in CLM) + age_days = clm%snowage / 86400.0_r8 + age_days = max(age_days, 0.001_r8) ! Prevent zero/negative values + + if (clm%t_grnd < tfrz) then + ! Cold/accumulating conditions - slow decay + albs = snal0 * (clm%albedo_accum_a ** (age_days ** 0.58_r8)) + albl = snal1 * (clm%albedo_accum_a ** (age_days ** 0.58_r8)) + else + ! Warm/melting conditions - fast decay + albs = snal0 * (clm%albedo_thaw_a ** (age_days ** 0.46_r8)) + albl = snal1 * (clm%albedo_thaw_a ** (age_days ** 0.46_r8)) + endif + + case (2) ! Tarboton - Arrhenius temperature dependence + ! ===================================================================== + ! Tarboton-style temperature-dependent albedo decay + ! Aging rate increases exponentially as temperature approaches freezing + ! Reference: Tarboton & Luce (1996), Utah Energy Balance Snow Model + ! ===================================================================== + + ! Temperature-dependent aging rate (faster near melting point) + ! Uses Arrhenius-type formulation + decay_factor = exp(5000.0_r8 * (1.0_r8/tfrz - 1.0_r8/max(clm%t_grnd, 200.0_r8))) + decay_factor = min(decay_factor, 10.0_r8) ! Cap to prevent extreme values + + ! Convert snow age to days + age_days = clm%snowage / 86400.0_r8 + + ! Modified age factor with temperature dependence + age = (age_days * decay_factor) / (1.0_r8 + age_days * decay_factor) + + albs = snal0 * (1.0_r8 - cons * age) + albl = snal1 * (1.0_r8 - conn * age) + + case default + ! Fallback to CLM default + age = 1._r8 - 1._r8/(1._r8 + clm%snowage) + albs = snal0*(1._r8 - cons*age) + albl = snal1*(1._r8 - conn*age) + + end select - age = 1.-1./(1.+clm%snowage) - albs = snal0*(1.-cons*age) - albl = snal1*(1.-conn*age) +! ========================================================================= +! Apply minimum albedo floor (prevents unrealistically low values) +! ========================================================================= + albs = max(clm%albedo_min, albs) + albl = max(clm%albedo_min, albl) - if (ind == 0) then +! ========================================================================= +! Zenith angle correction (applied to all schemes) +! ========================================================================= -! Czf corrects albedo of new snow for solar zenith + if (ind == 0) then - cff = ((1.+1./sl)/(1.+max(dble(0.001),coszen)*2.*sl )- 1./sl) - cff = max(cff,dble(0.)) - czf = 0.4*cff*(1.-albs) - albs = albs+czf - czf = 0.4*cff*(1.-albl) - albl = albl+czf + ! Czf corrects albedo of new snow for solar zenith + cff = ((1._r8 + 1._r8/sl)/(1._r8 + max(0.001_r8, coszen)*2._r8*sl) - 1._r8/sl) + cff = max(cff, 0.0_r8) + czf = 0.4_r8*cff*(1._r8 - albs) + albs = albs + czf + czf = 0.4_r8*cff*(1._r8 - albl) + albl = albl + czf endif diff --git a/pfsimulator/clm/clmtype.F90 b/pfsimulator/clm/clmtype.F90 index 7f61e6902..649627256 100644 --- a/pfsimulator/clm/clmtype.F90 +++ b/pfsimulator/clm/clmtype.F90 @@ -250,7 +250,44 @@ module clmtype real(r8) :: irr_threshold ! irrigation soil moisture threshold for deficit cycle @IMF integer :: threshold_type ! irrigation threshold type -- top layer, bottom layer, column avg. real(r8) :: irr_flag ! flag for irrigation or non-irrigation for a given day (based on threshold) - + +! Snow parameterization options @RMM 2025 + integer :: snow_partition_type ! rain-snow partition: 0=CLM, 1=wetbulb thresh, 2=wetbulb linear, 3=Dai, 4=Jennings + real(r8) :: tw_threshold ! wetbulb temperature threshold for snow [K], default 274.15 + real(r8) :: thin_snow_damping ! damping factor for thin snow energy [0-1], 0=off + real(r8) :: thin_snow_threshold ! SWE threshold for damping [kg/m2 or mm], default 50.0 + real(r8) :: snow_tcrit ! initial T classification threshold above tfrz [K], default 2.5 + real(r8) :: snow_t_low ! CLM method lower T threshold [K], default 273.16 + real(r8) :: snow_t_high ! CLM method upper T threshold [K], default 275.16 + real(r8) :: snow_transition_width ! WetbulbLinear half-width [K], default 1.0 + real(r8) :: dai_a ! Dai (2008) coefficient a, default -48.2292 + real(r8) :: dai_b ! Dai (2008) coefficient b, default 0.7205 + real(r8) :: dai_c ! Dai (2008) coefficient c, default 1.1662 + real(r8) :: dai_d ! Dai (2008) coefficient d, default 1.0223 + real(r8) :: jennings_a ! Jennings (2018) intercept, default -10.04 + real(r8) :: jennings_b ! Jennings (2018) T coefficient, default 1.41 + real(r8) :: jennings_g ! Jennings (2018) RH coefficient, default 0.09 + +! SZA-based snow melt damping parameters @RMM 2025 + real(r8) :: coszen ! current cosine of solar zenith angle (stored from clm_main) + real(r8) :: sza_snow_damping ! damping factor at high SZA [0-1], 1.0=disabled + real(r8) :: sza_damping_coszen_ref ! reference coszen below which damping applies, default 0.5 + real(r8) :: sza_damping_coszen_min ! coszen at which max damping applies, default 0.1 + +! Snow albedo parameterization options @RMM 2025 + integer :: albedo_scheme ! albedo scheme: 0=CLM, 1=VIC, 2=Tarboton + real(r8) :: albedo_vis_new ! fresh snow VIS albedo [0-1], default 0.95 + real(r8) :: albedo_nir_new ! fresh snow NIR albedo [0-1], default 0.65 + real(r8) :: albedo_min ! minimum albedo floor [0-1], default 0.4 + real(r8) :: albedo_decay_vis ! VIS decay coefficient [0-1], default 0.5 + real(r8) :: albedo_decay_nir ! NIR decay coefficient [0-1], default 0.2 + real(r8) :: albedo_accum_a ! VIC cold-phase decay base, default 0.94 + real(r8) :: albedo_thaw_a ! VIC melt-phase decay base, default 0.82 + +! Fractional snow covered area (frac_sno) options @RMM 2025 + integer :: frac_sno_type ! frac_sno scheme: 0=CLM (default), others TBD + real(r8) :: frac_sno_roughness ! roughness length for frac_sno [m], default=zlnd + real(r8) :: eflx_snomelt ! added to be consistent with lsm hybrid code real(r8) :: eflx_impsoil ! implicit evaporation for soil temperature equation (W/m**2) real(r8) :: eflx_lh_vege ! veg evaporation heat flux (W/m**2) [+ to atm] diff --git a/pfsimulator/clm/drv_getforce.F90 b/pfsimulator/clm/drv_getforce.F90 index ce9f4f477..f9b57cb89 100644 --- a/pfsimulator/clm/drv_getforce.F90 +++ b/pfsimulator/clm/drv_getforce.F90 @@ -33,7 +33,7 @@ subroutine drv_getforce (drv,tile,clm,nx,ny,sw_pf,lw_pf,prcp_pf,tas_pf,u_pf,v_pf use drv_module ! 1-D Land Model Driver variables use drv_tilemodule ! Tile-space variables use clmtype ! 1-D CLM variables - use clm_varcon, only : tfrz, tcrit + use clm_varcon, only : tfrz implicit none !=== Arguments =========================================================== @@ -112,8 +112,9 @@ subroutine drv_getforce (drv,tile,clm,nx,ny,sw_pf,lw_pf,prcp_pf,tas_pf,u_pf,v_pf !(Set upper limit of air temperature for snowfall at 275.65K. ! This cut-off was selected based on Fig. 1, Plate 3-1, of Snow ! Hydrology (1956)). + ! @RMM 2025: Now uses configurable snow_tcrit from clm1d struct if (prcp > 0.) then - if(clm(t)%forc_t > (tfrz + tcrit))then + if(clm(t)%forc_t > (tfrz + clm(t)%snow_tcrit))then clm(t)%itypprc = 1 clm(t)%forc_rain = prcp clm(t)%forc_snow = 0. diff --git a/pfsimulator/parflow_lib/parflow_proto_f.h b/pfsimulator/parflow_lib/parflow_proto_f.h index 61359969a..a57caeefb 100644 --- a/pfsimulator/parflow_lib/parflow_proto_f.h +++ b/pfsimulator/parflow_lib/parflow_proto_f.h @@ -120,31 +120,43 @@ void SADVECT(double *s, double *sn, #define CLM_LSM clm_lsm_ #endif -#define CALL_CLM_LSM(pressure_data, saturation_data, evap_trans_data, top, bottom, porosity_data, \ - dz_mult_data, istep, dt, t, start_time, dx, dy, dz, ix, iy, nx, ny, nz, \ - nx_f, ny_f, nz_f, nz_rz, ip, p, q, r, gnx, gny, rank, \ - sw_data, lw_data, prcp_data, tas_data, u_data, v_data, patm_data, qatm_data, \ - lai_data, sai_data, z0m_data, displa_data, \ - slope_x_data, slope_y_data, \ - eflx_lh_tot_data, eflx_lwrad_out_data, eflx_sh_tot_data, eflx_soil_grnd_data, \ - qflx_evap_tot_data, qflx_evap_grnd_data, qflx_evap_soi_data, qflx_evap_veg_data, qflx_tran_veg_data, \ - qflx_infl_data, swe_out_data, t_grnd_data, t_soil_data, \ - clm_dump_interval, clm_1d_out, clm_forc_veg, clm_file_dir, clm_file_dir_length, clm_bin_out_dir, write_CLM_binary, slope_accounting_CLM, \ - clm_beta_function, clm_veg_function, clm_veg_wilting, clm_veg_fieldc, clm_res_sat, \ - clm_irr_type, clm_irr_cycle, clm_irr_rate, clm_irr_start, clm_irr_stop, \ - clm_irr_threshold, qirr, qirr_inst, iflag, clm_irr_thresholdtype, soi_z, clm_next, clm_write_logs, clm_last_rst, clm_daily_rst, clm_water_stress_type, clm_nlevsoi, clm_nlevlak) \ - CLM_LSM(pressure_data, saturation_data, evap_trans_data, top, bottom, porosity_data, \ - dz_mult_data, &istep, &dt, &t, &start_time, &dx, &dy, &dz, &ix, &iy, &nx, &ny, &nz, &nx_f, &ny_f, &nz_f, &nz_rz, &ip, &p, &q, &r, &gnx, &gny, &rank, \ - sw_data, lw_data, prcp_data, tas_data, u_data, v_data, patm_data, qatm_data, \ - lai_data, sai_data, z0m_data, displa_data, \ - slope_x_data, slope_y_data, \ - eflx_lh_tot_data, eflx_lwrad_out_data, eflx_sh_tot_data, eflx_soil_grnd_data, \ - qflx_evap_tot_data, qflx_evap_grnd_data, qflx_evap_soi_data, qflx_evap_veg_data, qflx_tran_veg_data, \ - qflx_infl_data, swe_out_data, t_grnd_data, t_soil_data, \ - &clm_dump_interval, &clm_1d_out, &clm_forc_veg, clm_file_dir, &clm_file_dir_length, &clm_bin_out_dir, \ - &write_CLM_binary, &slope_accounting_CLM, &clm_beta_function, &clm_veg_function, &clm_veg_wilting, &clm_veg_fieldc, \ - &clm_res_sat, &clm_irr_type, &clm_irr_cycle, &clm_irr_rate, &clm_irr_start, &clm_irr_stop, \ - &clm_irr_threshold, qirr, qirr_inst, iflag, &clm_irr_thresholdtype, &soi_z, &clm_next, &clm_write_logs, &clm_last_rst, &clm_daily_rst, &clm_water_stress_type, &clm_nlevsoi, &clm_nlevlak); +#define CALL_CLM_LSM(pressure_data, saturation_data, evap_trans_data, top, bottom, porosity_data, \ + dz_mult_data, istep, dt, t, start_time, dx, dy, dz, ix, iy, nx, ny, nz, \ + nx_f, ny_f, nz_f, nz_rz, ip, p, q, r, gnx, gny, rank, \ + sw_data, lw_data, prcp_data, tas_data, u_data, v_data, patm_data, qatm_data, \ + lai_data, sai_data, z0m_data, displa_data, \ + slope_x_data, slope_y_data, \ + eflx_lh_tot_data, eflx_lwrad_out_data, eflx_sh_tot_data, eflx_soil_grnd_data, \ + qflx_evap_tot_data, qflx_evap_grnd_data, qflx_evap_soi_data, qflx_evap_veg_data, qflx_tran_veg_data, \ + qflx_infl_data, swe_out_data, t_grnd_data, t_soil_data, \ + clm_dump_interval, clm_1d_out, clm_forc_veg, clm_file_dir, clm_file_dir_length, clm_bin_out_dir, write_CLM_binary, slope_accounting_CLM, \ + clm_beta_function, clm_veg_function, clm_veg_wilting, clm_veg_fieldc, clm_res_sat, \ + clm_irr_type, clm_irr_cycle, clm_irr_rate, clm_irr_start, clm_irr_stop, \ + clm_irr_threshold, qirr, qirr_inst, iflag, clm_irr_thresholdtype, soi_z, clm_next, clm_write_logs, clm_last_rst, clm_daily_rst, clm_water_stress_type, clm_nlevsoi, clm_nlevlak, \ + clm_snow_partition, clm_tw_threshold, clm_thin_snow_damping, clm_thin_snow_threshold, \ + clm_snow_tcrit, clm_snow_t_low, clm_snow_t_high, clm_snow_transition_width, \ + clm_dai_a, clm_dai_b, clm_dai_c, clm_dai_d, clm_jennings_a, clm_jennings_b, clm_jennings_g, \ + clm_sza_snow_damping, clm_sza_damping_coszen_ref, clm_sza_damping_coszen_min, \ + clm_albedo_scheme, clm_albedo_vis_new, clm_albedo_nir_new, clm_albedo_min, clm_albedo_decay_vis, clm_albedo_decay_nir, clm_albedo_accum_a, clm_albedo_thaw_a, \ + clm_frac_sno_type, clm_frac_sno_roughness) \ + CLM_LSM(pressure_data, saturation_data, evap_trans_data, top, bottom, porosity_data, \ + dz_mult_data, &istep, &dt, &t, &start_time, &dx, &dy, &dz, &ix, &iy, &nx, &ny, &nz, &nx_f, &ny_f, &nz_f, &nz_rz, &ip, &p, &q, &r, &gnx, &gny, &rank, \ + sw_data, lw_data, prcp_data, tas_data, u_data, v_data, patm_data, qatm_data, \ + lai_data, sai_data, z0m_data, displa_data, \ + slope_x_data, slope_y_data, \ + eflx_lh_tot_data, eflx_lwrad_out_data, eflx_sh_tot_data, eflx_soil_grnd_data, \ + qflx_evap_tot_data, qflx_evap_grnd_data, qflx_evap_soi_data, qflx_evap_veg_data, qflx_tran_veg_data, \ + qflx_infl_data, swe_out_data, t_grnd_data, t_soil_data, \ + &clm_dump_interval, &clm_1d_out, &clm_forc_veg, clm_file_dir, &clm_file_dir_length, &clm_bin_out_dir, \ + &write_CLM_binary, &slope_accounting_CLM, &clm_beta_function, &clm_veg_function, &clm_veg_wilting, &clm_veg_fieldc, \ + &clm_res_sat, &clm_irr_type, &clm_irr_cycle, &clm_irr_rate, &clm_irr_start, &clm_irr_stop, \ + &clm_irr_threshold, qirr, qirr_inst, iflag, &clm_irr_thresholdtype, &soi_z, &clm_next, &clm_write_logs, &clm_last_rst, &clm_daily_rst, &clm_water_stress_type, &clm_nlevsoi, &clm_nlevlak, \ + &clm_snow_partition, &clm_tw_threshold, &clm_thin_snow_damping, &clm_thin_snow_threshold, \ + &clm_snow_tcrit, &clm_snow_t_low, &clm_snow_t_high, &clm_snow_transition_width, \ + &clm_dai_a, &clm_dai_b, &clm_dai_c, &clm_dai_d, &clm_jennings_a, &clm_jennings_b, &clm_jennings_g, \ + &clm_sza_snow_damping, &clm_sza_damping_coszen_ref, &clm_sza_damping_coszen_min, \ + &clm_albedo_scheme, &clm_albedo_vis_new, &clm_albedo_nir_new, &clm_albedo_min, &clm_albedo_decay_vis, &clm_albedo_decay_nir, &clm_albedo_accum_a, &clm_albedo_thaw_a, \ + &clm_frac_sno_type, &clm_frac_sno_roughness); void CLM_LSM(double *pressure_data, double *saturation_data, double *evap_trans_data, double *top, double *bottom, double *porosity_data, double *dz_mult_data, int *istep, double *dt, double *t, double *start_time, @@ -160,7 +172,15 @@ void CLM_LSM(double *pressure_data, double *saturation_data, double *evap_trans_ int *clm_veg_function, double *clm_veg_wilting, double *clm_veg_fieldc, double *clm_res_sat, int *clm_irr_type, int *clm_irr_cycle, double *clm_irr_rate, double *clm_irr_start, double *clm_irr_stop, double *clm_irr_threshold, double *qirr, double *qirr_inst, double *iflag, int *clm_irr_thresholdtype, int *soi_z, - int *clm_next, int *clm_write_logs, int *clm_last_rst, int *clm_daily_rst, int *clm_water_stress_type, int *clm_nlevsoi, int *clm_nlevlak); + int *clm_next, int *clm_write_logs, int *clm_last_rst, int *clm_daily_rst, int *clm_water_stress_type, int *clm_nlevsoi, int *clm_nlevlak, + int *clm_snow_partition, double *clm_tw_threshold, double *clm_thin_snow_damping, double *clm_thin_snow_threshold, + double *clm_snow_tcrit, double *clm_snow_t_low, double *clm_snow_t_high, double *clm_snow_transition_width, + double *clm_dai_a, double *clm_dai_b, double *clm_dai_c, double *clm_dai_d, + double *clm_jennings_a, double *clm_jennings_b, double *clm_jennings_g, + double *clm_sza_snow_damping, double *clm_sza_damping_coszen_ref, double *clm_sza_damping_coszen_min, + int *clm_albedo_scheme, double *clm_albedo_vis_new, double *clm_albedo_nir_new, double *clm_albedo_min, + double *clm_albedo_decay_vis, double *clm_albedo_decay_nir, double *clm_albedo_accum_a, double *clm_albedo_thaw_a, + int *clm_frac_sno_type, double *clm_frac_sno_roughness); /* @RMM CRUNCHFLOW.F90*/ //#define CRUNCHFLOW crunchflow_ diff --git a/pfsimulator/parflow_lib/solver_richards.c b/pfsimulator/parflow_lib/solver_richards.c index 2b9d5c4d1..1add8c20d 100644 --- a/pfsimulator/parflow_lib/solver_richards.c +++ b/pfsimulator/parflow_lib/solver_richards.c @@ -197,6 +197,42 @@ typedef struct { double clm_irr_threshold; /* CLM irrigation schedule -- soil moisture threshold for deficit cycle */ int clm_irr_thresholdtype; /* Deficit-based saturation criteria (top, bottom, column avg) */ + /* Snow parameterization options @RMM 2025 */ + int clm_snow_partition; /* CLM snow partition type: 0=CLM, 1=wetbulb thresh, 2=wetbulb linear, 3=Dai, 4=Jennings */ + double clm_tw_threshold; /* CLM wetbulb temperature threshold for snow [K] */ + double clm_thin_snow_damping; /* CLM thin snow energy damping factor [0-1] */ + double clm_thin_snow_threshold; /* CLM SWE threshold for damping [kg/m2] */ + double clm_snow_tcrit; /* Initial T classification threshold above tfrz [K], default 2.5 */ + double clm_snow_t_low; /* CLM method lower T threshold [K], default 273.16 */ + double clm_snow_t_high; /* CLM method upper T threshold [K], default 275.16 */ + double clm_snow_transition_width; /* WetbulbLinear half-width [K], default 1.0 */ + double clm_dai_a; /* Dai (2008) coefficient a, default -48.2292 */ + double clm_dai_b; /* Dai (2008) coefficient b, default 0.7205 */ + double clm_dai_c; /* Dai (2008) coefficient c, default 1.1662 */ + double clm_dai_d; /* Dai (2008) coefficient d, default 1.0223 */ + double clm_jennings_a; /* Jennings (2018) intercept, default -10.04 */ + double clm_jennings_b; /* Jennings (2018) T coefficient, default 1.41 */ + double clm_jennings_g; /* Jennings (2018) RH coefficient, default 0.09 */ + + /* SZA-based snow damping parameters @RMM 2025 */ + double clm_sza_snow_damping; /* SZA damping factor [0-1], 1.0=disabled */ + double clm_sza_damping_coszen_ref; /* Reference coszen for damping onset (0.5 = 60 deg) */ + double clm_sza_damping_coszen_min; /* Coszen at maximum damping (0.1 = 84 deg) */ + + /* Snow albedo parameterization options @RMM 2025 */ + int clm_albedo_scheme; /* CLM albedo scheme: 0=CLM, 1=VIC, 2=Tarboton */ + double clm_albedo_vis_new; /* Fresh snow VIS albedo [0-1] */ + double clm_albedo_nir_new; /* Fresh snow NIR albedo [0-1] */ + double clm_albedo_min; /* Minimum albedo floor [0-1] */ + double clm_albedo_decay_vis; /* VIS decay coefficient [0-1] */ + double clm_albedo_decay_nir; /* NIR decay coefficient [0-1] */ + double clm_albedo_accum_a; /* VIC cold-phase decay base */ + double clm_albedo_thaw_a; /* VIC melt-phase decay base */ + + /* Fractional snow covered area (frac_sno) options @RMM 2025 */ + int clm_frac_sno_type; /* frac_sno scheme: 0=CLM (default), others TBD */ + double clm_frac_sno_roughness; /* roughness length for frac_sno [m], default=0.01 (zlnd) */ + int clm_reuse_count; /* NBE: Number of times to use each CLM input */ int clm_write_logs; /* NBE: Write the processor logs for CLM or not */ int clm_last_rst; /* NBE: Only write/overwrite one rst file or write a lot of them */ @@ -993,7 +1029,8 @@ SetupRichards(PFModule * this_module) if (public_xtra->write_silo_qx_overland || public_xtra->print_qx_overland || public_xtra->write_pdi_qx_overland - || public_xtra->write_netcdf_qx_overland) + || public_xtra->write_netcdf_qx_overland + || public_xtra->surface_predictor) { instance_xtra->q_overlnd_x = NewVectorType(grid2d, 1, 1, vector_cell_centered_2D); @@ -1007,7 +1044,8 @@ SetupRichards(PFModule * this_module) if (public_xtra->write_silo_qy_overland || public_xtra->print_qy_overland || public_xtra->write_pdi_qy_overland - || public_xtra->write_netcdf_qy_overland) + || public_xtra->write_netcdf_qy_overland + || public_xtra->surface_predictor) { instance_xtra->q_overlnd_y = NewVectorType(grid2d, 1, 1, vector_cell_centered_2D); @@ -2765,7 +2803,35 @@ AdvanceRichards(PFModule * this_module, double start_time, /* Starting time clm_daily_rst, clm_water_stress_type, public_xtra->clm_nz, - public_xtra->clm_nz); + public_xtra->clm_nz, + public_xtra->clm_snow_partition, + public_xtra->clm_tw_threshold, + public_xtra->clm_thin_snow_damping, + public_xtra->clm_thin_snow_threshold, + public_xtra->clm_snow_tcrit, + public_xtra->clm_snow_t_low, + public_xtra->clm_snow_t_high, + public_xtra->clm_snow_transition_width, + public_xtra->clm_dai_a, + public_xtra->clm_dai_b, + public_xtra->clm_dai_c, + public_xtra->clm_dai_d, + public_xtra->clm_jennings_a, + public_xtra->clm_jennings_b, + public_xtra->clm_jennings_g, + public_xtra->clm_sza_snow_damping, + public_xtra->clm_sza_damping_coszen_ref, + public_xtra->clm_sza_damping_coszen_min, + public_xtra->clm_albedo_scheme, + public_xtra->clm_albedo_vis_new, + public_xtra->clm_albedo_nir_new, + public_xtra->clm_albedo_min, + public_xtra->clm_albedo_decay_vis, + public_xtra->clm_albedo_decay_nir, + public_xtra->clm_albedo_accum_a, + public_xtra->clm_albedo_thaw_a, + public_xtra->clm_frac_sno_type, + public_xtra->clm_frac_sno_roughness); break; } @@ -5599,6 +5665,181 @@ SolverRichardsNewPublicXtra(char *name) sprintf(key, "%s.CLM.FieldCapacity", name); public_xtra->clm_veg_fieldc = GetDoubleDefault(key, 1.0); + /* @RMM 2025 Snow parameterization options */ + NameArray snow_switch_na; + snow_switch_na = NA_NewNameArray("CLM WetbulbThreshold WetbulbLinear Dai Jennings"); + sprintf(key, "%s.CLM.SnowPartition", name); + switch_name = GetStringDefault(key, "CLM"); + switch_value = NA_NameToIndexExitOnError(snow_switch_na, switch_name, key); + switch (switch_value) + { + case 0: + { + public_xtra->clm_snow_partition = 0; + break; + } + + case 1: + { + public_xtra->clm_snow_partition = 1; + break; + } + + case 2: + { + public_xtra->clm_snow_partition = 2; + break; + } + + case 3: + { + public_xtra->clm_snow_partition = 3; + break; + } + + case 4: + { + public_xtra->clm_snow_partition = 4; + break; + } + + default: + { + InputError("Invalid switch value <%s> for key <%s>", switch_name, key); + } + } + NA_FreeNameArray(snow_switch_na); + + sprintf(key, "%s.CLM.WetbulbThreshold", name); + public_xtra->clm_tw_threshold = GetDoubleDefault(key, 274.15); + + sprintf(key, "%s.CLM.SnowTCrit", name); + public_xtra->clm_snow_tcrit = GetDoubleDefault(key, 2.5); + + sprintf(key, "%s.CLM.SnowTLow", name); + public_xtra->clm_snow_t_low = GetDoubleDefault(key, 273.16); + + sprintf(key, "%s.CLM.SnowTHigh", name); + public_xtra->clm_snow_t_high = GetDoubleDefault(key, 275.16); + + sprintf(key, "%s.CLM.SnowTransitionWidth", name); + public_xtra->clm_snow_transition_width = GetDoubleDefault(key, 1.0); + + sprintf(key, "%s.CLM.DaiCoeffA", name); + public_xtra->clm_dai_a = GetDoubleDefault(key, -48.2292); + + sprintf(key, "%s.CLM.DaiCoeffB", name); + public_xtra->clm_dai_b = GetDoubleDefault(key, 0.7205); + + sprintf(key, "%s.CLM.DaiCoeffC", name); + public_xtra->clm_dai_c = GetDoubleDefault(key, 1.1662); + + sprintf(key, "%s.CLM.DaiCoeffD", name); + public_xtra->clm_dai_d = GetDoubleDefault(key, 1.0223); + + sprintf(key, "%s.CLM.JenningsCoeffA", name); + public_xtra->clm_jennings_a = GetDoubleDefault(key, -10.04); + + sprintf(key, "%s.CLM.JenningsCoeffB", name); + public_xtra->clm_jennings_b = GetDoubleDefault(key, 1.41); + + sprintf(key, "%s.CLM.JenningsCoeffG", name); + public_xtra->clm_jennings_g = GetDoubleDefault(key, 0.09); + + sprintf(key, "%s.CLM.ThinSnowDamping", name); + public_xtra->clm_thin_snow_damping = GetDoubleDefault(key, 0.0); + + sprintf(key, "%s.CLM.ThinSnowThreshold", name); + public_xtra->clm_thin_snow_threshold = GetDoubleDefault(key, 50.0); + + /* @RMM 2025 SZA-based snow melt damping */ + sprintf(key, "%s.CLM.SZASnowDamping", name); + public_xtra->clm_sza_snow_damping = GetDoubleDefault(key, 1.0); + + sprintf(key, "%s.CLM.SZADampingCoszenRef", name); + public_xtra->clm_sza_damping_coszen_ref = GetDoubleDefault(key, 0.5); + + sprintf(key, "%s.CLM.SZADampingCoszenMin", name); + public_xtra->clm_sza_damping_coszen_min = GetDoubleDefault(key, 0.1); + + /* @RMM 2025 Snow albedo parameterization options */ + NameArray albedo_switch_na; + albedo_switch_na = NA_NewNameArray("CLM VIC Tarboton"); + sprintf(key, "%s.CLM.AlbedoScheme", name); + switch_name = GetStringDefault(key, "CLM"); + switch_value = NA_NameToIndexExitOnError(albedo_switch_na, switch_name, key); + switch (switch_value) + { + case 0: + { + public_xtra->clm_albedo_scheme = 0; + break; + } + + case 1: + { + public_xtra->clm_albedo_scheme = 1; + break; + } + + case 2: + { + public_xtra->clm_albedo_scheme = 2; + break; + } + + default: + { + InputError("Invalid switch value <%s> for key <%s>", switch_name, key); + } + } + NA_FreeNameArray(albedo_switch_na); + + sprintf(key, "%s.CLM.AlbedoVisNew", name); + public_xtra->clm_albedo_vis_new = GetDoubleDefault(key, 0.95); + + sprintf(key, "%s.CLM.AlbedoNirNew", name); + public_xtra->clm_albedo_nir_new = GetDoubleDefault(key, 0.65); + + sprintf(key, "%s.CLM.AlbedoMin", name); + public_xtra->clm_albedo_min = GetDoubleDefault(key, 0.4); + + sprintf(key, "%s.CLM.AlbedoDecayVis", name); + public_xtra->clm_albedo_decay_vis = GetDoubleDefault(key, 0.5); + + sprintf(key, "%s.CLM.AlbedoDecayNir", name); + public_xtra->clm_albedo_decay_nir = GetDoubleDefault(key, 0.2); + + sprintf(key, "%s.CLM.AlbedoAccumA", name); + public_xtra->clm_albedo_accum_a = GetDoubleDefault(key, 0.94); + + sprintf(key, "%s.CLM.AlbedoThawA", name); + public_xtra->clm_albedo_thaw_a = GetDoubleDefault(key, 0.82); + + /* @RMM 2025 Fractional snow covered area (frac_sno) options */ + NameArray frac_sno_switch_na; + frac_sno_switch_na = NA_NewNameArray("CLM"); + sprintf(key, "%s.CLM.FracSnoScheme", name); + switch_name = GetStringDefault(key, "CLM"); + switch_value = NA_NameToIndexExitOnError(frac_sno_switch_na, switch_name, key); + switch (switch_value) + { + case 0: + { + public_xtra->clm_frac_sno_type = 0; + break; + } + + default: + { + InputError("Invalid switch value <%s> for key <%s>", switch_name, key); + } + } + NA_FreeNameArray(frac_sno_switch_na); + + sprintf(key, "%s.CLM.FracSnoRoughness", name); + public_xtra->clm_frac_sno_roughness = GetDoubleDefault(key, 0.01); + /* IMF Write CLM as Silo (default=False) */ sprintf(key, "%s.WriteSiloCLM", name); switch_name = GetStringDefault(key, "False"); diff --git a/pftools/python/parflow/tools/io.py b/pftools/python/parflow/tools/io.py index 387925a72..882b71fef 100644 --- a/pftools/python/parflow/tools/io.py +++ b/pftools/python/parflow/tools/io.py @@ -1491,7 +1491,7 @@ def _read_vegm(file_name): in the vegm.dat file except for x/y """ # Assume first two lines are comments and use generic column names - df = pd.read_csv(file_name, delim_whitespace=True, skiprows=2, header=None) + df = pd.read_csv(file_name, sep="\s+", skiprows=2, header=None) df.columns = [f"c{i}" for i in range(df.shape[1])] # Number of columns and rows determined by last line of file diff --git a/pftools/python/pyproject.toml b/pftools/python/pyproject.toml index 5a807443f..b50bcffc2 100644 --- a/pftools/python/pyproject.toml +++ b/pftools/python/pyproject.toml @@ -34,11 +34,11 @@ all = [ "imageio>=2.9.0", "h5py", "twine", - "black" + "black==25.12.0" ] pfsol = ["imageio>=2.9.0"] pdi = ["h5py"] -dev = ["twine", "black"] +dev = ["twine", "black==25.12.0"] [project.urls] Homepage = "https://github.com/parflow/parflow/tree/master/pftools/python" diff --git a/pftools/python/requirements_dev.txt b/pftools/python/requirements_dev.txt index 46c7b2ef5..525404890 100644 --- a/pftools/python/requirements_dev.txt +++ b/pftools/python/requirements_dev.txt @@ -1,3 +1,3 @@ ###### Requirements for development of PFTools package ###### twine -black +black==25.12.0 diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.alpha.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.alpha.pfb new file mode 100644 index 000000000..4d9bbe3b6 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.alpha.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.clm_output.00001.C.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.clm_output.00001.C.pfb new file mode 100644 index 000000000..913ab7f61 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.clm_output.00001.C.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.clm_output.00002.C.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.clm_output.00002.C.pfb new file mode 100644 index 000000000..a9b8df1c7 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.clm_output.00002.C.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.clm_output.00003.C.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.clm_output.00003.C.pfb new file mode 100644 index 000000000..8408cfce5 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.clm_output.00003.C.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.clm_output.00004.C.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.clm_output.00004.C.pfb new file mode 100644 index 000000000..e6be9cce2 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.clm_output.00004.C.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.clm_output.00005.C.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.clm_output.00005.C.pfb new file mode 100644 index 000000000..1f5243261 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.clm_output.00005.C.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.mask.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.mask.pfb new file mode 100644 index 000000000..caa1ffcc9 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.mask.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.n.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.n.pfb new file mode 100644 index 000000000..18472f609 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.n.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.perm_x.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.perm_x.pfb new file mode 100644 index 000000000..614ae8457 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.perm_x.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.perm_y.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.perm_y.pfb new file mode 100644 index 000000000..614ae8457 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.perm_y.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.perm_z.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.perm_z.pfb new file mode 100644 index 000000000..614ae8457 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.perm_z.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.porosity.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.porosity.pfb new file mode 100644 index 000000000..13c2a6121 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.porosity.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.press.00000.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.press.00000.pfb new file mode 100644 index 000000000..fb17a1ccd Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.press.00000.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.press.00001.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.press.00001.pfb new file mode 100644 index 000000000..5a12e0176 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.press.00001.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.press.00002.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.press.00002.pfb new file mode 100644 index 000000000..bbd842e8b Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.press.00002.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.press.00003.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.press.00003.pfb new file mode 100644 index 000000000..285dc361e Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.press.00003.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.press.00004.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.press.00004.pfb new file mode 100644 index 000000000..5df72d5e6 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.press.00004.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.press.00005.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.press.00005.pfb new file mode 100644 index 000000000..9aef45d14 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.press.00005.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.satur.00000.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.satur.00000.pfb new file mode 100644 index 000000000..8ea84826c Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.satur.00000.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.satur.00001.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.satur.00001.pfb new file mode 100644 index 000000000..e586e4604 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.satur.00001.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.satur.00002.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.satur.00002.pfb new file mode 100644 index 000000000..fd9005443 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.satur.00002.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.satur.00003.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.satur.00003.pfb new file mode 100644 index 000000000..93c3dbc8e Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.satur.00003.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.satur.00004.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.satur.00004.pfb new file mode 100644 index 000000000..59ede8041 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.satur.00004.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.satur.00005.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.satur.00005.pfb new file mode 100644 index 000000000..94f7e95c6 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.satur.00005.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.specific_storage.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.specific_storage.pfb new file mode 100644 index 000000000..1297279b9 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.specific_storage.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.sres.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.sres.pfb new file mode 100644 index 000000000..5175c79ba Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.sres.pfb differ diff --git a/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.ssat.pfb b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.ssat.pfb new file mode 100644 index 000000000..caa1ffcc9 Binary files /dev/null and b/test/correct_output/clm_snow_albedo/clm_snow_albedo.out.ssat.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.alpha.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.alpha.pfb new file mode 100644 index 000000000..4d9bbe3b6 Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.alpha.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.clm_output.00001.C.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.clm_output.00001.C.pfb new file mode 100644 index 000000000..94798175e Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.clm_output.00001.C.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.clm_output.00002.C.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.clm_output.00002.C.pfb new file mode 100644 index 000000000..5c93d3de8 Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.clm_output.00002.C.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.clm_output.00003.C.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.clm_output.00003.C.pfb new file mode 100644 index 000000000..6a6c90913 Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.clm_output.00003.C.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.clm_output.00004.C.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.clm_output.00004.C.pfb new file mode 100644 index 000000000..0f0b9ce3f Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.clm_output.00004.C.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.clm_output.00005.C.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.clm_output.00005.C.pfb new file mode 100644 index 000000000..e3d26590c Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.clm_output.00005.C.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.mask.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.mask.pfb new file mode 100644 index 000000000..caa1ffcc9 Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.mask.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.n.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.n.pfb new file mode 100644 index 000000000..18472f609 Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.n.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.perm_x.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.perm_x.pfb new file mode 100644 index 000000000..614ae8457 Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.perm_x.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.perm_y.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.perm_y.pfb new file mode 100644 index 000000000..614ae8457 Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.perm_y.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.perm_z.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.perm_z.pfb new file mode 100644 index 000000000..614ae8457 Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.perm_z.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.porosity.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.porosity.pfb new file mode 100644 index 000000000..13c2a6121 Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.porosity.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.press.00000.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.press.00000.pfb new file mode 100644 index 000000000..fb17a1ccd Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.press.00000.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.press.00001.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.press.00001.pfb new file mode 100644 index 000000000..a4b4243ad Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.press.00001.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.press.00002.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.press.00002.pfb new file mode 100644 index 000000000..2efd5efb0 Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.press.00002.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.press.00003.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.press.00003.pfb new file mode 100644 index 000000000..823c25d4c Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.press.00003.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.press.00004.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.press.00004.pfb new file mode 100644 index 000000000..364e0935e Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.press.00004.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.press.00005.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.press.00005.pfb new file mode 100644 index 000000000..56700bd9a Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.press.00005.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.satur.00000.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.satur.00000.pfb new file mode 100644 index 000000000..8ea84826c Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.satur.00000.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.satur.00001.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.satur.00001.pfb new file mode 100644 index 000000000..df1930a62 Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.satur.00001.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.satur.00002.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.satur.00002.pfb new file mode 100644 index 000000000..8bb2139c0 Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.satur.00002.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.satur.00003.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.satur.00003.pfb new file mode 100644 index 000000000..aef234fcc Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.satur.00003.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.satur.00004.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.satur.00004.pfb new file mode 100644 index 000000000..d190a7194 Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.satur.00004.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.satur.00005.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.satur.00005.pfb new file mode 100644 index 000000000..1c2f7abc2 Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.satur.00005.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.specific_storage.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.specific_storage.pfb new file mode 100644 index 000000000..1297279b9 Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.specific_storage.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.sres.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.sres.pfb new file mode 100644 index 000000000..5175c79ba Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.sres.pfb differ diff --git a/test/correct_output/clm_snow_dai/clm_snow_dai.out.ssat.pfb b/test/correct_output/clm_snow_dai/clm_snow_dai.out.ssat.pfb new file mode 100644 index 000000000..caa1ffcc9 Binary files /dev/null and b/test/correct_output/clm_snow_dai/clm_snow_dai.out.ssat.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.alpha.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.alpha.pfb new file mode 100644 index 000000000..4d9bbe3b6 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.alpha.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.clm_output.00001.C.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.clm_output.00001.C.pfb new file mode 100644 index 000000000..913ab7f61 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.clm_output.00001.C.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.clm_output.00002.C.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.clm_output.00002.C.pfb new file mode 100644 index 000000000..a9b8df1c7 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.clm_output.00002.C.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.clm_output.00003.C.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.clm_output.00003.C.pfb new file mode 100644 index 000000000..8408cfce5 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.clm_output.00003.C.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.clm_output.00004.C.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.clm_output.00004.C.pfb new file mode 100644 index 000000000..e6be9cce2 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.clm_output.00004.C.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.clm_output.00005.C.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.clm_output.00005.C.pfb new file mode 100644 index 000000000..1f5243261 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.clm_output.00005.C.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.mask.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.mask.pfb new file mode 100644 index 000000000..caa1ffcc9 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.mask.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.n.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.n.pfb new file mode 100644 index 000000000..18472f609 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.n.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.perm_x.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.perm_x.pfb new file mode 100644 index 000000000..614ae8457 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.perm_x.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.perm_y.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.perm_y.pfb new file mode 100644 index 000000000..614ae8457 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.perm_y.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.perm_z.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.perm_z.pfb new file mode 100644 index 000000000..614ae8457 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.perm_z.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.porosity.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.porosity.pfb new file mode 100644 index 000000000..13c2a6121 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.porosity.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.press.00000.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.press.00000.pfb new file mode 100644 index 000000000..fb17a1ccd Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.press.00000.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.press.00001.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.press.00001.pfb new file mode 100644 index 000000000..5a12e0176 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.press.00001.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.press.00002.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.press.00002.pfb new file mode 100644 index 000000000..bbd842e8b Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.press.00002.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.press.00003.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.press.00003.pfb new file mode 100644 index 000000000..285dc361e Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.press.00003.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.press.00004.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.press.00004.pfb new file mode 100644 index 000000000..5df72d5e6 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.press.00004.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.press.00005.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.press.00005.pfb new file mode 100644 index 000000000..9aef45d14 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.press.00005.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.satur.00000.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.satur.00000.pfb new file mode 100644 index 000000000..8ea84826c Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.satur.00000.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.satur.00001.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.satur.00001.pfb new file mode 100644 index 000000000..e586e4604 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.satur.00001.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.satur.00002.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.satur.00002.pfb new file mode 100644 index 000000000..fd9005443 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.satur.00002.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.satur.00003.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.satur.00003.pfb new file mode 100644 index 000000000..93c3dbc8e Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.satur.00003.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.satur.00004.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.satur.00004.pfb new file mode 100644 index 000000000..59ede8041 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.satur.00004.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.satur.00005.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.satur.00005.pfb new file mode 100644 index 000000000..94f7e95c6 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.satur.00005.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.specific_storage.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.specific_storage.pfb new file mode 100644 index 000000000..1297279b9 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.specific_storage.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.sres.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.sres.pfb new file mode 100644 index 000000000..5175c79ba Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.sres.pfb differ diff --git a/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.ssat.pfb b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.ssat.pfb new file mode 100644 index 000000000..caa1ffcc9 Binary files /dev/null and b/test/correct_output/clm_snow_defaults/clm_snow_defaults.out.ssat.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.alpha.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.alpha.pfb new file mode 100644 index 000000000..4d9bbe3b6 Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.alpha.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.clm_output.00001.C.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.clm_output.00001.C.pfb new file mode 100644 index 000000000..18a379ee4 Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.clm_output.00001.C.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.clm_output.00002.C.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.clm_output.00002.C.pfb new file mode 100644 index 000000000..0a08c15df Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.clm_output.00002.C.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.clm_output.00003.C.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.clm_output.00003.C.pfb new file mode 100644 index 000000000..0b0d09e20 Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.clm_output.00003.C.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.clm_output.00004.C.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.clm_output.00004.C.pfb new file mode 100644 index 000000000..6da55cabc Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.clm_output.00004.C.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.clm_output.00005.C.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.clm_output.00005.C.pfb new file mode 100644 index 000000000..e30dd556b Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.clm_output.00005.C.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.mask.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.mask.pfb new file mode 100644 index 000000000..caa1ffcc9 Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.mask.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.n.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.n.pfb new file mode 100644 index 000000000..18472f609 Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.n.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.perm_x.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.perm_x.pfb new file mode 100644 index 000000000..614ae8457 Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.perm_x.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.perm_y.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.perm_y.pfb new file mode 100644 index 000000000..614ae8457 Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.perm_y.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.perm_z.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.perm_z.pfb new file mode 100644 index 000000000..614ae8457 Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.perm_z.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.porosity.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.porosity.pfb new file mode 100644 index 000000000..13c2a6121 Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.porosity.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.press.00000.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.press.00000.pfb new file mode 100644 index 000000000..fb17a1ccd Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.press.00000.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.press.00001.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.press.00001.pfb new file mode 100644 index 000000000..1a9b10938 Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.press.00001.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.press.00002.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.press.00002.pfb new file mode 100644 index 000000000..1dd4eb886 Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.press.00002.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.press.00003.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.press.00003.pfb new file mode 100644 index 000000000..0f734e50d Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.press.00003.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.press.00004.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.press.00004.pfb new file mode 100644 index 000000000..680210e01 Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.press.00004.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.press.00005.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.press.00005.pfb new file mode 100644 index 000000000..c075db80a Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.press.00005.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.satur.00000.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.satur.00000.pfb new file mode 100644 index 000000000..8ea84826c Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.satur.00000.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.satur.00001.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.satur.00001.pfb new file mode 100644 index 000000000..33bea6a2d Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.satur.00001.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.satur.00002.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.satur.00002.pfb new file mode 100644 index 000000000..698810f18 Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.satur.00002.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.satur.00003.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.satur.00003.pfb new file mode 100644 index 000000000..f8a475223 Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.satur.00003.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.satur.00004.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.satur.00004.pfb new file mode 100644 index 000000000..f1378c6ab Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.satur.00004.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.satur.00005.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.satur.00005.pfb new file mode 100644 index 000000000..8ddb7318f Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.satur.00005.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.specific_storage.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.specific_storage.pfb new file mode 100644 index 000000000..1297279b9 Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.specific_storage.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.sres.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.sres.pfb new file mode 100644 index 000000000..5175c79ba Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.sres.pfb differ diff --git a/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.ssat.pfb b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.ssat.pfb new file mode 100644 index 000000000..caa1ffcc9 Binary files /dev/null and b/test/correct_output/clm_snow_jennings/clm_snow_jennings.out.ssat.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.alpha.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.alpha.pfb new file mode 100644 index 000000000..4d9bbe3b6 Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.alpha.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.clm_output.00001.C.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.clm_output.00001.C.pfb new file mode 100644 index 000000000..6fc86420b Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.clm_output.00001.C.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.clm_output.00002.C.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.clm_output.00002.C.pfb new file mode 100644 index 000000000..aef4ac729 Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.clm_output.00002.C.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.clm_output.00003.C.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.clm_output.00003.C.pfb new file mode 100644 index 000000000..5cacbb91b Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.clm_output.00003.C.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.clm_output.00004.C.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.clm_output.00004.C.pfb new file mode 100644 index 000000000..f6a1b2420 Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.clm_output.00004.C.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.clm_output.00005.C.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.clm_output.00005.C.pfb new file mode 100644 index 000000000..aed1a0a91 Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.clm_output.00005.C.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.mask.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.mask.pfb new file mode 100644 index 000000000..caa1ffcc9 Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.mask.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.n.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.n.pfb new file mode 100644 index 000000000..18472f609 Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.n.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.perm_x.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.perm_x.pfb new file mode 100644 index 000000000..614ae8457 Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.perm_x.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.perm_y.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.perm_y.pfb new file mode 100644 index 000000000..614ae8457 Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.perm_y.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.perm_z.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.perm_z.pfb new file mode 100644 index 000000000..614ae8457 Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.perm_z.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.porosity.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.porosity.pfb new file mode 100644 index 000000000..13c2a6121 Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.porosity.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.press.00000.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.press.00000.pfb new file mode 100644 index 000000000..fb17a1ccd Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.press.00000.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.press.00001.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.press.00001.pfb new file mode 100644 index 000000000..3d49e69b6 Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.press.00001.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.press.00002.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.press.00002.pfb new file mode 100644 index 000000000..bcc3694c4 Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.press.00002.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.press.00003.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.press.00003.pfb new file mode 100644 index 000000000..1a64c6294 Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.press.00003.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.press.00004.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.press.00004.pfb new file mode 100644 index 000000000..82db23ef5 Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.press.00004.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.press.00005.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.press.00005.pfb new file mode 100644 index 000000000..b2c81266f Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.press.00005.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.satur.00000.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.satur.00000.pfb new file mode 100644 index 000000000..8ea84826c Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.satur.00000.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.satur.00001.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.satur.00001.pfb new file mode 100644 index 000000000..8d0b55d86 Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.satur.00001.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.satur.00002.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.satur.00002.pfb new file mode 100644 index 000000000..9a199137e Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.satur.00002.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.satur.00003.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.satur.00003.pfb new file mode 100644 index 000000000..167dc3220 Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.satur.00003.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.satur.00004.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.satur.00004.pfb new file mode 100644 index 000000000..4fbd18571 Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.satur.00004.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.satur.00005.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.satur.00005.pfb new file mode 100644 index 000000000..0f62cb2ce Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.satur.00005.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.specific_storage.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.specific_storage.pfb new file mode 100644 index 000000000..1297279b9 Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.specific_storage.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.sres.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.sres.pfb new file mode 100644 index 000000000..5175c79ba Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.sres.pfb differ diff --git a/test/correct_output/clm_snow_partition/clm_snow_partition.out.ssat.pfb b/test/correct_output/clm_snow_partition/clm_snow_partition.out.ssat.pfb new file mode 100644 index 000000000..caa1ffcc9 Binary files /dev/null and b/test/correct_output/clm_snow_partition/clm_snow_partition.out.ssat.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.alpha.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.alpha.pfb new file mode 100644 index 000000000..4d9bbe3b6 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.alpha.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.clm_output.00001.C.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.clm_output.00001.C.pfb new file mode 100644 index 000000000..0e311ace8 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.clm_output.00001.C.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.clm_output.00002.C.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.clm_output.00002.C.pfb new file mode 100644 index 000000000..2d1237435 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.clm_output.00002.C.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.clm_output.00003.C.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.clm_output.00003.C.pfb new file mode 100644 index 000000000..3e8b97b12 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.clm_output.00003.C.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.clm_output.00004.C.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.clm_output.00004.C.pfb new file mode 100644 index 000000000..2e86d5008 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.clm_output.00004.C.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.clm_output.00005.C.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.clm_output.00005.C.pfb new file mode 100644 index 000000000..c01978cc5 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.clm_output.00005.C.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.mask.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.mask.pfb new file mode 100644 index 000000000..caa1ffcc9 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.mask.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.n.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.n.pfb new file mode 100644 index 000000000..18472f609 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.n.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.perm_x.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.perm_x.pfb new file mode 100644 index 000000000..614ae8457 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.perm_x.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.perm_y.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.perm_y.pfb new file mode 100644 index 000000000..614ae8457 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.perm_y.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.perm_z.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.perm_z.pfb new file mode 100644 index 000000000..614ae8457 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.perm_z.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.porosity.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.porosity.pfb new file mode 100644 index 000000000..13c2a6121 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.porosity.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.press.00000.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.press.00000.pfb new file mode 100644 index 000000000..fb17a1ccd Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.press.00000.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.press.00001.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.press.00001.pfb new file mode 100644 index 000000000..6f7a24ef0 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.press.00001.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.press.00002.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.press.00002.pfb new file mode 100644 index 000000000..d7be441b4 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.press.00002.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.press.00003.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.press.00003.pfb new file mode 100644 index 000000000..2ffbade56 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.press.00003.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.press.00004.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.press.00004.pfb new file mode 100644 index 000000000..8968a8d98 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.press.00004.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.press.00005.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.press.00005.pfb new file mode 100644 index 000000000..8f4890c29 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.press.00005.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.satur.00000.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.satur.00000.pfb new file mode 100644 index 000000000..8ea84826c Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.satur.00000.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.satur.00001.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.satur.00001.pfb new file mode 100644 index 000000000..ad5d27fc8 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.satur.00001.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.satur.00002.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.satur.00002.pfb new file mode 100644 index 000000000..9a50c030c Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.satur.00002.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.satur.00003.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.satur.00003.pfb new file mode 100644 index 000000000..352a75773 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.satur.00003.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.satur.00004.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.satur.00004.pfb new file mode 100644 index 000000000..39e44c62b Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.satur.00004.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.satur.00005.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.satur.00005.pfb new file mode 100644 index 000000000..94732cfcc Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.satur.00005.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.specific_storage.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.specific_storage.pfb new file mode 100644 index 000000000..1297279b9 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.specific_storage.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.sres.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.sres.pfb new file mode 100644 index 000000000..5175c79ba Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.sres.pfb differ diff --git a/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.ssat.pfb b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.ssat.pfb new file mode 100644 index 000000000..caa1ffcc9 Binary files /dev/null and b/test/correct_output/clm_snow_sza_damping/clm_snow_sza_damping.out.ssat.pfb differ diff --git a/test/python/clm/CMakeLists.txt b/test/python/clm/CMakeLists.txt index b57c71bb8..2ed70d6e2 100644 --- a/test/python/clm/CMakeLists.txt +++ b/test/python/clm/CMakeLists.txt @@ -7,7 +7,13 @@ if (${PARFLOW_HAVE_CLM}) if (${PARFLOW_HAVE_HYPRE}) list(APPEND CLM_TESTS clm-reuse - clm_rz_water_stress) + clm_rz_water_stress + clm_snow_defaults + clm_snow_partition + clm_snow_albedo + clm_snow_dai + clm_snow_jennings + clm_snow_sza_damping) list(APPEND CLM_2D_TESTS clm_slope) diff --git a/test/python/clm/clm_snow_albedo.py b/test/python/clm/clm_snow_albedo.py new file mode 100644 index 000000000..9028d6f6e --- /dev/null +++ b/test/python/clm/clm_snow_albedo.py @@ -0,0 +1,413 @@ +# ----------------------------------------------------------------------------- +# CLM Snow Albedo Test +# +# This test verifies that the snow albedo scheme options work correctly. +# It runs with the Tarboton albedo scheme and compares to reference output. +# +# The Tarboton scheme uses Arrhenius temperature-dependent aging where +# decay accelerates near the melting point. +# +# Reference: Tarboton & Luce (1996) WRR, doi:10.1029/96WR01049 +# ----------------------------------------------------------------------------- + +import sys +import argparse + +from parflow import Run +from parflow.tools.fs import mkdir, cp, get_absolute_path, rm +from parflow.tools.compare import pf_test_file + +run_name = "clm_snow_albedo" +clm = Run(run_name, __file__) + +# ----------------------------------------------------------------------------- +# Making output directories and copying input files +# ----------------------------------------------------------------------------- + +dir_name = get_absolute_path("test_output/" + run_name) +mkdir(dir_name) + +directories = [ + "qflx_evap_grnd", + "eflx_lh_tot", + "qflx_evap_tot", + "qflx_tran_veg", + "correct_output", + "qflx_infl", + "swe_out", + "eflx_lwrad_out", + "t_grnd", + "diag_out", + "qflx_evap_soi", + "eflx_soil_grnd", + "eflx_sh_tot", + "qflx_evap_veg", + "qflx_top_soil", +] + +for directory in directories: + mkdir(dir_name + "/" + directory) + +cp("$PF_SRC/test/tcl/clm/drv_clmin.dat", dir_name) +cp("$PF_SRC/test/tcl/clm/drv_vegm.dat", dir_name) +cp("$PF_SRC/test/tcl/clm/drv_vegp.dat", dir_name) +cp("$PF_SRC/test/tcl/clm/snow_forcing.1hr.txt", dir_name) + +# ----------------------------------------------------------------------------- +# File input version number +# ----------------------------------------------------------------------------- + +clm.FileVersion = 4 + +# ----------------------------------------------------------------------------- +# Process Topology +# ----------------------------------------------------------------------------- + +parser = argparse.ArgumentParser() +parser.add_argument("-p", "--p", default=1) +parser.add_argument("-q", "--q", default=1) +parser.add_argument("-r", "--r", default=1) +args = parser.parse_args() + +clm.Process.Topology.P = args.p +clm.Process.Topology.Q = args.q +clm.Process.Topology.R = args.r + +# ----------------------------------------------------------------------------- +# Computational Grid +# ----------------------------------------------------------------------------- + +clm.ComputationalGrid.Lower.X = 0.0 +clm.ComputationalGrid.Lower.Y = 0.0 +clm.ComputationalGrid.Lower.Z = 0.0 + +clm.ComputationalGrid.DX = 1000.0 +clm.ComputationalGrid.DY = 1000.0 +clm.ComputationalGrid.DZ = 0.5 + +clm.ComputationalGrid.NX = 5 +clm.ComputationalGrid.NY = 5 +clm.ComputationalGrid.NZ = 10 + +# ----------------------------------------------------------------------------- +# The Names of the GeomInputs +# ----------------------------------------------------------------------------- + +clm.GeomInput.Names = "domain_input" + +# ----------------------------------------------------------------------------- +# Domain Geometry Input +# ----------------------------------------------------------------------------- + +clm.GeomInput.domain_input.InputType = "Box" +clm.GeomInput.domain_input.GeomName = "domain" + +# ----------------------------------------------------------------------------- +# Domain Geometry +# ----------------------------------------------------------------------------- + +clm.Geom.domain.Lower.X = 0.0 +clm.Geom.domain.Lower.Y = 0.0 +clm.Geom.domain.Lower.Z = 0.0 + +clm.Geom.domain.Upper.X = 5000.0 +clm.Geom.domain.Upper.Y = 5000.0 +clm.Geom.domain.Upper.Z = 5.0 + +clm.Geom.domain.Patches = "x_lower x_upper y_lower y_upper z_lower z_upper" + +# ----------------------------------------------------------------------------- +# Perm +# ----------------------------------------------------------------------------- + +clm.Geom.Perm.Names = "domain" +clm.Geom.domain.Perm.Type = "Constant" +clm.Geom.domain.Perm.Value = 0.2 + +clm.Perm.TensorType = "TensorByGeom" +clm.Geom.Perm.TensorByGeom.Names = "domain" + +clm.Geom.domain.Perm.TensorValX = 1.0 +clm.Geom.domain.Perm.TensorValY = 1.0 +clm.Geom.domain.Perm.TensorValZ = 1.0 + +# ----------------------------------------------------------------------------- +# Specific Storage +# ----------------------------------------------------------------------------- + +clm.SpecificStorage.Type = "Constant" +clm.SpecificStorage.GeomNames = "domain" +clm.Geom.domain.SpecificStorage.Value = 1.0e-6 + +# ----------------------------------------------------------------------------- +# Phases +# ----------------------------------------------------------------------------- + +clm.Phase.Names = "water" +clm.Phase.water.Density.Type = "Constant" +clm.Phase.water.Density.Value = 1.0 +clm.Phase.water.Viscosity.Type = "Constant" +clm.Phase.water.Viscosity.Value = 1.0 + +# ----------------------------------------------------------------------------- +# Contaminants +# ----------------------------------------------------------------------------- + +clm.Contaminants.Names = "" + +# ----------------------------------------------------------------------------- +# Gravity +# ----------------------------------------------------------------------------- + +clm.Gravity = 1.0 + +# ----------------------------------------------------------------------------- +# Setup timing info +# ----------------------------------------------------------------------------- + +clm.TimingInfo.BaseUnit = 1.0 +clm.TimingInfo.StartCount = 0 +clm.TimingInfo.StartTime = 0.0 +clm.TimingInfo.StopTime = 5 +clm.TimingInfo.DumpInterval = -1 +clm.TimeStep.Type = "Constant" +clm.TimeStep.Value = 1.0 + +# ----------------------------------------------------------------------------- +# Porosity +# ----------------------------------------------------------------------------- + +clm.Geom.Porosity.GeomNames = "domain" +clm.Geom.domain.Porosity.Type = "Constant" +clm.Geom.domain.Porosity.Value = 0.390 + +# ----------------------------------------------------------------------------- +# Domain +# ----------------------------------------------------------------------------- + +clm.Domain.GeomName = "domain" + +# ----------------------------------------------------------------------------- +# Mobility +# ----------------------------------------------------------------------------- + +clm.Phase.water.Mobility.Type = "Constant" +clm.Phase.water.Mobility.Value = 1.0 + +# ----------------------------------------------------------------------------- +# Relative Permeability +# ----------------------------------------------------------------------------- + +clm.Phase.RelPerm.Type = "VanGenuchten" +clm.Phase.RelPerm.GeomNames = "domain" +clm.Geom.domain.RelPerm.Alpha = 3.5 +clm.Geom.domain.RelPerm.N = 2.0 + +# ----------------------------------------------------------------------------- +# Saturation +# ----------------------------------------------------------------------------- + +clm.Phase.Saturation.Type = "VanGenuchten" +clm.Phase.Saturation.GeomNames = "domain" +clm.Geom.domain.Saturation.Alpha = 3.5 +clm.Geom.domain.Saturation.N = 2.0 +clm.Geom.domain.Saturation.SRes = 0.01 +clm.Geom.domain.Saturation.SSat = 1.0 + +# ----------------------------------------------------------------------------- +# Wells +# ----------------------------------------------------------------------------- + +clm.Wells.Names = "" + +# ----------------------------------------------------------------------------- +# Time Cycles +# ----------------------------------------------------------------------------- + +clm.Cycle.Names = "constant" +clm.Cycle.constant.Names = "alltime" +clm.Cycle.constant.alltime.Length = 1 +clm.Cycle.constant.Repeat = -1 + +# ----------------------------------------------------------------------------- +# Boundary Conditions: Pressure +# ----------------------------------------------------------------------------- + +clm.BCPressure.PatchNames = "x_lower x_upper y_lower y_upper z_lower z_upper" + +clm.Patch.x_lower.BCPressure.Type = "FluxConst" +clm.Patch.x_lower.BCPressure.Cycle = "constant" +clm.Patch.x_lower.BCPressure.alltime.Value = 0.0 + +clm.Patch.y_lower.BCPressure.Type = "FluxConst" +clm.Patch.y_lower.BCPressure.Cycle = "constant" +clm.Patch.y_lower.BCPressure.alltime.Value = 0.0 + +clm.Patch.z_lower.BCPressure.Type = "FluxConst" +clm.Patch.z_lower.BCPressure.Cycle = "constant" +clm.Patch.z_lower.BCPressure.alltime.Value = 0.0 + +clm.Patch.x_upper.BCPressure.Type = "FluxConst" +clm.Patch.x_upper.BCPressure.Cycle = "constant" +clm.Patch.x_upper.BCPressure.alltime.Value = 0.0 + +clm.Patch.y_upper.BCPressure.Type = "FluxConst" +clm.Patch.y_upper.BCPressure.Cycle = "constant" +clm.Patch.y_upper.BCPressure.alltime.Value = 0.0 + +clm.Patch.z_upper.BCPressure.Type = "OverlandFlow" +clm.Patch.z_upper.BCPressure.Cycle = "constant" +clm.Patch.z_upper.BCPressure.alltime.Value = 0.0 + +# ----------------------------------------------------------------------------- +# Topo slopes in x-direction +# ----------------------------------------------------------------------------- + +clm.TopoSlopesX.Type = "Constant" +clm.TopoSlopesX.GeomNames = "domain" +clm.TopoSlopesX.Geom.domain.Value = -0.001 + +# ----------------------------------------------------------------------------- +# Topo slopes in y-direction +# ----------------------------------------------------------------------------- + +clm.TopoSlopesY.Type = "Constant" +clm.TopoSlopesY.GeomNames = "domain" +clm.TopoSlopesY.Geom.domain.Value = 0.001 + +# ----------------------------------------------------------------------------- +# Mannings coefficient +# ----------------------------------------------------------------------------- + +clm.Mannings.Type = "Constant" +clm.Mannings.GeomNames = "domain" +clm.Mannings.Geom.domain.Value = 5.52e-6 + +# ----------------------------------------------------------------------------- +# Phase sources +# ----------------------------------------------------------------------------- + +clm.PhaseSources.water.Type = "Constant" +clm.PhaseSources.water.GeomNames = "domain" +clm.PhaseSources.water.Geom.domain.Value = 0.0 + +# ----------------------------------------------------------------------------- +# Exact solution specification for error calculations +# ----------------------------------------------------------------------------- + +clm.KnownSolution = "NoKnownSolution" + +# ----------------------------------------------------------------------------- +# Set solver parameters +# ----------------------------------------------------------------------------- + +clm.Solver = "Richards" +clm.Solver.MaxIter = 500 + +clm.Solver.Nonlinear.MaxIter = 15 +clm.Solver.Nonlinear.ResidualTol = 1e-9 +clm.Solver.Nonlinear.EtaChoice = "EtaConstant" +clm.Solver.Nonlinear.EtaValue = 0.01 +clm.Solver.Nonlinear.UseJacobian = True +clm.Solver.Nonlinear.StepTol = 1e-20 +clm.Solver.Nonlinear.Globalization = "LineSearch" +clm.Solver.Linear.KrylovDimension = 15 +clm.Solver.Linear.MaxRestart = 2 + +clm.Solver.Linear.Preconditioner = "PFMG" +clm.Solver.PrintSubsurf = False +clm.Solver.Drop = 1e-20 +clm.Solver.AbsTol = 1e-9 + +# ----------------------------------------------------------------------------- +# CLM Settings +# ----------------------------------------------------------------------------- + +clm.Solver.LSM = "CLM" +clm.Solver.CLM.MetForcing = "1D" +clm.Solver.CLM.MetFileName = "snow_forcing.1hr.txt" +clm.Solver.CLM.MetFilePath = "." + +clm.Solver.WriteSiloCLM = False +clm.Solver.WriteSiloEvapTrans = False +clm.Solver.WriteSiloOverlandBCFlux = False +clm.Solver.PrintCLM = True + +clm.Solver.CLM.Print1dOut = False +clm.Solver.BinaryOutDir = False +clm.Solver.WriteCLMBinary = False +clm.Solver.CLM.CLMDumpInterval = 1 +clm.Solver.CLM.WriteLogs = False +clm.Solver.CLM.WriteLastRST = True +clm.Solver.CLM.DailyRST = True +clm.Solver.CLM.SingleFile = True + +# ----------------------------------------------------------------------------- +# Snow Albedo Scheme - Use Tarboton temperature-dependent aging +# This exercises the Tarboton albedo code path +# ----------------------------------------------------------------------------- + +clm.Solver.CLM.AlbedoScheme = "Tarboton" +clm.Solver.CLM.AlbedoVisNew = 0.95 +clm.Solver.CLM.AlbedoNirNew = 0.65 +clm.Solver.CLM.AlbedoMin = 0.4 +clm.Solver.CLM.AlbedoDecayVis = 0.5 +clm.Solver.CLM.AlbedoDecayNir = 0.2 + +# ----------------------------------------------------------------------------- +# Initial conditions: water pressure +# ----------------------------------------------------------------------------- + +clm.ICPressure.Type = "HydroStaticPatch" +clm.ICPressure.GeomNames = "domain" +clm.Geom.domain.ICPressure.Value = -2.0 +clm.Geom.domain.ICPressure.RefGeom = "domain" +clm.Geom.domain.ICPressure.RefPatch = "z_upper" + +# ----------------------------------------------------------------------------- +# Run and Unload the ParFlow output files +# ----------------------------------------------------------------------------- + +correct_output_dir_name = get_absolute_path("../../correct_output/" + run_name) + +clm.run(working_directory=dir_name) + +# ----------------------------------------------------------------------------- +# Tests - Compare pressure, saturation, and CLM output +# ----------------------------------------------------------------------------- + +passed = True + +for i in range(6): + timestep = str(i).rjust(5, "0") + filename = f"/{run_name}.out.press.{timestep}.pfb" + if not pf_test_file( + dir_name + filename, + correct_output_dir_name + filename, + f"Max difference in Pressure for timestep {timestep}", + ): + passed = False + filename = f"/{run_name}.out.satur.{timestep}.pfb" + if not pf_test_file( + dir_name + filename, + correct_output_dir_name + filename, + f"Max difference in Saturation for timestep {timestep}", + ): + passed = False + + if i > 0: + filename = f"/{run_name}.out.clm_output.{timestep}.C.pfb" + if not pf_test_file( + dir_name + filename, + correct_output_dir_name + filename, + f"Max difference in CLM output for timestep {timestep}", + ): + passed = False + +rm(dir_name) + +if passed: + print(f"{run_name} : PASSED") +else: + print(f"{run_name} : FAILED") + sys.exit(1) diff --git a/test/python/clm/clm_snow_dai.py b/test/python/clm/clm_snow_dai.py new file mode 100644 index 000000000..39fae060e --- /dev/null +++ b/test/python/clm/clm_snow_dai.py @@ -0,0 +1,425 @@ +# ----------------------------------------------------------------------------- +# CLM Snow Dai Test +# +# This test verifies that the Dai (2008) rain-snow partitioning method works. +# Formula: F(%) = a * [tanh(b*(T-c)) - d] +# +# The Dai method uses a hyperbolic tangent function fitted to 30 years of +# global weather station observations of precipitation phase. +# +# Reference: Dai (2008) GRL, doi:10.1029/2008GL033295 +# +# New keys tested (defaults from Table 1a, Land, ANN): +# Solver.CLM.SnowPartition = "Dai" +# Solver.CLM.DaiCoeffA = -48.2292 (scaling factor) +# Solver.CLM.DaiCoeffB = 0.7205 (slope) +# Solver.CLM.DaiCoeffC = 1.1662 (half-frequency T in C) +# Solver.CLM.DaiCoeffD = 1.0223 (asymmetry) +# ----------------------------------------------------------------------------- + +import sys +import argparse + +from parflow import Run +from parflow.tools.fs import mkdir, cp, get_absolute_path, rm +from parflow.tools.compare import pf_test_file + +run_name = "clm_snow_dai" +clm = Run(run_name, __file__) + +# ----------------------------------------------------------------------------- +# Making output directories and copying input files +# ----------------------------------------------------------------------------- + +dir_name = get_absolute_path("test_output/" + run_name) +mkdir(dir_name) + +directories = [ + "qflx_evap_grnd", + "eflx_lh_tot", + "qflx_evap_tot", + "qflx_tran_veg", + "correct_output", + "qflx_infl", + "swe_out", + "eflx_lwrad_out", + "t_grnd", + "diag_out", + "qflx_evap_soi", + "eflx_soil_grnd", + "eflx_sh_tot", + "qflx_evap_veg", + "qflx_top_soil", +] + +for directory in directories: + mkdir(dir_name + "/" + directory) + +cp("$PF_SRC/test/tcl/clm/drv_clmin.dat", dir_name) +cp("$PF_SRC/test/tcl/clm/drv_vegm.dat", dir_name) +cp("$PF_SRC/test/tcl/clm/drv_vegp.dat", dir_name) +cp("$PF_SRC/test/tcl/clm/snow_forcing.1hr.txt", dir_name) + +# ----------------------------------------------------------------------------- +# File input version number +# ----------------------------------------------------------------------------- + +clm.FileVersion = 4 + +# ----------------------------------------------------------------------------- +# Process Topology +# ----------------------------------------------------------------------------- + +parser = argparse.ArgumentParser() +parser.add_argument("-p", "--p", default=1) +parser.add_argument("-q", "--q", default=1) +parser.add_argument("-r", "--r", default=1) +args = parser.parse_args() + +clm.Process.Topology.P = args.p +clm.Process.Topology.Q = args.q +clm.Process.Topology.R = args.r + +# ----------------------------------------------------------------------------- +# Computational Grid +# ----------------------------------------------------------------------------- + +clm.ComputationalGrid.Lower.X = 0.0 +clm.ComputationalGrid.Lower.Y = 0.0 +clm.ComputationalGrid.Lower.Z = 0.0 + +clm.ComputationalGrid.DX = 1000.0 +clm.ComputationalGrid.DY = 1000.0 +clm.ComputationalGrid.DZ = 0.5 + +clm.ComputationalGrid.NX = 5 +clm.ComputationalGrid.NY = 5 +clm.ComputationalGrid.NZ = 10 + +# ----------------------------------------------------------------------------- +# The Names of the GeomInputs +# ----------------------------------------------------------------------------- + +clm.GeomInput.Names = "domain_input" + +# ----------------------------------------------------------------------------- +# Domain Geometry Input +# ----------------------------------------------------------------------------- + +clm.GeomInput.domain_input.InputType = "Box" +clm.GeomInput.domain_input.GeomName = "domain" + +# ----------------------------------------------------------------------------- +# Domain Geometry +# ----------------------------------------------------------------------------- + +clm.Geom.domain.Lower.X = 0.0 +clm.Geom.domain.Lower.Y = 0.0 +clm.Geom.domain.Lower.Z = 0.0 + +clm.Geom.domain.Upper.X = 5000.0 +clm.Geom.domain.Upper.Y = 5000.0 +clm.Geom.domain.Upper.Z = 5.0 + +clm.Geom.domain.Patches = "x_lower x_upper y_lower y_upper z_lower z_upper" + +# ----------------------------------------------------------------------------- +# Perm +# ----------------------------------------------------------------------------- + +clm.Geom.Perm.Names = "domain" +clm.Geom.domain.Perm.Type = "Constant" +clm.Geom.domain.Perm.Value = 0.2 + +clm.Perm.TensorType = "TensorByGeom" +clm.Geom.Perm.TensorByGeom.Names = "domain" + +clm.Geom.domain.Perm.TensorValX = 1.0 +clm.Geom.domain.Perm.TensorValY = 1.0 +clm.Geom.domain.Perm.TensorValZ = 1.0 + +# ----------------------------------------------------------------------------- +# Specific Storage +# ----------------------------------------------------------------------------- + +clm.SpecificStorage.Type = "Constant" +clm.SpecificStorage.GeomNames = "domain" +clm.Geom.domain.SpecificStorage.Value = 1.0e-6 + +# ----------------------------------------------------------------------------- +# Phases +# ----------------------------------------------------------------------------- + +clm.Phase.Names = "water" +clm.Phase.water.Density.Type = "Constant" +clm.Phase.water.Density.Value = 1.0 +clm.Phase.water.Viscosity.Type = "Constant" +clm.Phase.water.Viscosity.Value = 1.0 + +# ----------------------------------------------------------------------------- +# Contaminants +# ----------------------------------------------------------------------------- + +clm.Contaminants.Names = "" + +# ----------------------------------------------------------------------------- +# Gravity +# ----------------------------------------------------------------------------- + +clm.Gravity = 1.0 + +# ----------------------------------------------------------------------------- +# Setup timing info +# ----------------------------------------------------------------------------- + +clm.TimingInfo.BaseUnit = 1.0 +clm.TimingInfo.StartCount = 0 +clm.TimingInfo.StartTime = 0.0 +clm.TimingInfo.StopTime = 5 +clm.TimingInfo.DumpInterval = -1 +clm.TimeStep.Type = "Constant" +clm.TimeStep.Value = 1.0 + +# ----------------------------------------------------------------------------- +# Porosity +# ----------------------------------------------------------------------------- + +clm.Geom.Porosity.GeomNames = "domain" +clm.Geom.domain.Porosity.Type = "Constant" +clm.Geom.domain.Porosity.Value = 0.390 + +# ----------------------------------------------------------------------------- +# Domain +# ----------------------------------------------------------------------------- + +clm.Domain.GeomName = "domain" + +# ----------------------------------------------------------------------------- +# Mobility +# ----------------------------------------------------------------------------- + +clm.Phase.water.Mobility.Type = "Constant" +clm.Phase.water.Mobility.Value = 1.0 + +# ----------------------------------------------------------------------------- +# Relative Permeability +# ----------------------------------------------------------------------------- + +clm.Phase.RelPerm.Type = "VanGenuchten" +clm.Phase.RelPerm.GeomNames = "domain" +clm.Geom.domain.RelPerm.Alpha = 3.5 +clm.Geom.domain.RelPerm.N = 2.0 + +# ----------------------------------------------------------------------------- +# Saturation +# ----------------------------------------------------------------------------- + +clm.Phase.Saturation.Type = "VanGenuchten" +clm.Phase.Saturation.GeomNames = "domain" +clm.Geom.domain.Saturation.Alpha = 3.5 +clm.Geom.domain.Saturation.N = 2.0 +clm.Geom.domain.Saturation.SRes = 0.01 +clm.Geom.domain.Saturation.SSat = 1.0 + +# ----------------------------------------------------------------------------- +# Wells +# ----------------------------------------------------------------------------- + +clm.Wells.Names = "" + +# ----------------------------------------------------------------------------- +# Time Cycles +# ----------------------------------------------------------------------------- + +clm.Cycle.Names = "constant" +clm.Cycle.constant.Names = "alltime" +clm.Cycle.constant.alltime.Length = 1 +clm.Cycle.constant.Repeat = -1 + +# ----------------------------------------------------------------------------- +# Boundary Conditions: Pressure +# ----------------------------------------------------------------------------- + +clm.BCPressure.PatchNames = "x_lower x_upper y_lower y_upper z_lower z_upper" + +clm.Patch.x_lower.BCPressure.Type = "FluxConst" +clm.Patch.x_lower.BCPressure.Cycle = "constant" +clm.Patch.x_lower.BCPressure.alltime.Value = 0.0 + +clm.Patch.y_lower.BCPressure.Type = "FluxConst" +clm.Patch.y_lower.BCPressure.Cycle = "constant" +clm.Patch.y_lower.BCPressure.alltime.Value = 0.0 + +clm.Patch.z_lower.BCPressure.Type = "FluxConst" +clm.Patch.z_lower.BCPressure.Cycle = "constant" +clm.Patch.z_lower.BCPressure.alltime.Value = 0.0 + +clm.Patch.x_upper.BCPressure.Type = "FluxConst" +clm.Patch.x_upper.BCPressure.Cycle = "constant" +clm.Patch.x_upper.BCPressure.alltime.Value = 0.0 + +clm.Patch.y_upper.BCPressure.Type = "FluxConst" +clm.Patch.y_upper.BCPressure.Cycle = "constant" +clm.Patch.y_upper.BCPressure.alltime.Value = 0.0 + +clm.Patch.z_upper.BCPressure.Type = "OverlandFlow" +clm.Patch.z_upper.BCPressure.Cycle = "constant" +clm.Patch.z_upper.BCPressure.alltime.Value = 0.0 + +# ----------------------------------------------------------------------------- +# Topo slopes in x-direction +# ----------------------------------------------------------------------------- + +clm.TopoSlopesX.Type = "Constant" +clm.TopoSlopesX.GeomNames = "domain" +clm.TopoSlopesX.Geom.domain.Value = -0.001 + +# ----------------------------------------------------------------------------- +# Topo slopes in y-direction +# ----------------------------------------------------------------------------- + +clm.TopoSlopesY.Type = "Constant" +clm.TopoSlopesY.GeomNames = "domain" +clm.TopoSlopesY.Geom.domain.Value = 0.001 + +# ----------------------------------------------------------------------------- +# Mannings coefficient +# ----------------------------------------------------------------------------- + +clm.Mannings.Type = "Constant" +clm.Mannings.GeomNames = "domain" +clm.Mannings.Geom.domain.Value = 5.52e-6 + +# ----------------------------------------------------------------------------- +# Phase sources +# ----------------------------------------------------------------------------- + +clm.PhaseSources.water.Type = "Constant" +clm.PhaseSources.water.GeomNames = "domain" +clm.PhaseSources.water.Geom.domain.Value = 0.0 + +# ----------------------------------------------------------------------------- +# Exact solution specification for error calculations +# ----------------------------------------------------------------------------- + +clm.KnownSolution = "NoKnownSolution" + +# ----------------------------------------------------------------------------- +# Set solver parameters +# ----------------------------------------------------------------------------- + +clm.Solver = "Richards" +clm.Solver.MaxIter = 500 + +clm.Solver.Nonlinear.MaxIter = 15 +clm.Solver.Nonlinear.ResidualTol = 1e-9 +clm.Solver.Nonlinear.EtaChoice = "EtaConstant" +clm.Solver.Nonlinear.EtaValue = 0.01 +clm.Solver.Nonlinear.UseJacobian = True +clm.Solver.Nonlinear.StepTol = 1e-20 +clm.Solver.Nonlinear.Globalization = "LineSearch" +clm.Solver.Linear.KrylovDimension = 15 +clm.Solver.Linear.MaxRestart = 2 + +clm.Solver.Linear.Preconditioner = "PFMG" +clm.Solver.PrintSubsurf = False +clm.Solver.Drop = 1e-20 +clm.Solver.AbsTol = 1e-9 + +# ----------------------------------------------------------------------------- +# CLM Settings +# ----------------------------------------------------------------------------- + +clm.Solver.LSM = "CLM" +clm.Solver.CLM.MetForcing = "1D" +clm.Solver.CLM.MetFileName = "snow_forcing.1hr.txt" +clm.Solver.CLM.MetFilePath = "." + +clm.Solver.WriteSiloCLM = False +clm.Solver.WriteSiloEvapTrans = False +clm.Solver.WriteSiloOverlandBCFlux = False +clm.Solver.PrintCLM = True + +clm.Solver.CLM.Print1dOut = False +clm.Solver.BinaryOutDir = False +clm.Solver.WriteCLMBinary = False +clm.Solver.CLM.CLMDumpInterval = 1 +clm.Solver.CLM.WriteLogs = False +clm.Solver.CLM.WriteLastRST = True +clm.Solver.CLM.DailyRST = True +clm.Solver.CLM.SingleFile = True + +# ----------------------------------------------------------------------------- +# Snow Partitioning - Use Dai (2008) sigmoidal method +# F(%) = a * [tanh(b*(T-c)) - d], converted to fraction internally +# Coefficients from Dai (2008) Table 1a (Land, ANN): +# a = -48.2292 (scaling factor, negative) +# b = 0.7205 (slope at half-frequency temperature) +# c = 1.1662 (half-frequency temperature in C, ~50% snow) +# d = 1.0223 (asymmetry parameter) +# Reference: Dai (2008) GRL doi:10.1029/2008GL033295 +# ----------------------------------------------------------------------------- + +clm.Solver.CLM.SnowPartition = "Dai" +clm.Solver.CLM.DaiCoeffA = -48.2292 +clm.Solver.CLM.DaiCoeffB = 0.7205 +clm.Solver.CLM.DaiCoeffC = 1.1662 +clm.Solver.CLM.DaiCoeffD = 1.0223 + +# ----------------------------------------------------------------------------- +# Initial conditions: water pressure +# ----------------------------------------------------------------------------- + +clm.ICPressure.Type = "HydroStaticPatch" +clm.ICPressure.GeomNames = "domain" +clm.Geom.domain.ICPressure.Value = -2.0 +clm.Geom.domain.ICPressure.RefGeom = "domain" +clm.Geom.domain.ICPressure.RefPatch = "z_upper" + +# ----------------------------------------------------------------------------- +# Run ParFlow +# ----------------------------------------------------------------------------- + +clm.run(working_directory=dir_name) + +# ----------------------------------------------------------------------------- +# Tests - Compare pressure, saturation, and CLM output +# ----------------------------------------------------------------------------- + +correct_output_dir_name = get_absolute_path("$PF_SRC/test/correct_output/" + run_name) + +passed = True + +for i in range(6): + timestep = str(i).rjust(5, "0") + filename = f"/{run_name}.out.press.{timestep}.pfb" + if not pf_test_file( + dir_name + filename, + correct_output_dir_name + filename, + f"Max difference in Pressure for timestep {timestep}", + ): + passed = False + filename = f"/{run_name}.out.satur.{timestep}.pfb" + if not pf_test_file( + dir_name + filename, + correct_output_dir_name + filename, + f"Max difference in Saturation for timestep {timestep}", + ): + passed = False + + if i > 0: + filename = f"/{run_name}.out.clm_output.{timestep}.C.pfb" + if not pf_test_file( + dir_name + filename, + correct_output_dir_name + filename, + f"Max difference in CLM output for timestep {timestep}", + ): + passed = False + +rm(dir_name) + +if passed: + print(f"{run_name} : PASSED") +else: + print(f"{run_name} : FAILED") + sys.exit(1) diff --git a/test/python/clm/clm_snow_defaults.py b/test/python/clm/clm_snow_defaults.py new file mode 100644 index 000000000..8e8d977e6 --- /dev/null +++ b/test/python/clm/clm_snow_defaults.py @@ -0,0 +1,444 @@ +# ----------------------------------------------------------------------------- +# CLM Snow Defaults - Backward Compatibility Test +# +# This test verifies that when all new CLM snow parameterization keys +# are left at their default values, output is identical to unmodified CLM. +# This ensures backward compatibility for existing workflows. +# +# New keys tested (all at defaults): +# Rain-snow partitioning: +# Solver.CLM.SnowPartition = "CLM" +# Solver.CLM.WetbulbThreshold = 274.15 +# Solver.CLM.SnowTCrit = 2.5 +# Solver.CLM.SnowTLow = 273.16 +# Solver.CLM.SnowTHigh = 275.16 +# Solver.CLM.SnowTransitionWidth = 1.0 +# Solver.CLM.DaiCoeffA/B/C/D (Dai 2008 method) +# Solver.CLM.JenningsCoeffA/B/G (Jennings 2018 method) +# Melt damping: +# Solver.CLM.ThinSnowDamping = 0.0 +# Solver.CLM.ThinSnowThreshold = 50.0 +# Solver.CLM.SZASnowDamping = 1.0 (disabled) +# Solver.CLM.SZADampingCoszenRef = 0.5 +# Solver.CLM.SZADampingCoszenMin = 0.1 +# Albedo: +# Solver.CLM.AlbedoScheme = "CLM" +# Solver.CLM.AlbedoVisNew = 0.95 +# Solver.CLM.AlbedoNirNew = 0.65 +# Solver.CLM.AlbedoMin = 0.4 +# Solver.CLM.AlbedoDecayVis = 0.5 +# Solver.CLM.AlbedoDecayNir = 0.2 +# Solver.CLM.AlbedoAccumA = 0.94 +# Solver.CLM.AlbedoThawA = 0.82 +# Fractional snow cover: +# Solver.CLM.FracSnoScheme = "CLM" +# Solver.CLM.FracSnoRoughness = 0.01 +# ----------------------------------------------------------------------------- + +import sys +import argparse + +from parflow import Run +from parflow.tools.fs import mkdir, cp, get_absolute_path, rm +from parflow.tools.compare import pf_test_file + +run_name = "clm_snow_defaults" +clm = Run(run_name, __file__) + +# ----------------------------------------------------------------------------- +# Making output directories and copying input files +# ----------------------------------------------------------------------------- + +dir_name = get_absolute_path("test_output/" + run_name) +mkdir(dir_name) + +directories = [ + "qflx_evap_grnd", + "eflx_lh_tot", + "qflx_evap_tot", + "qflx_tran_veg", + "correct_output", + "qflx_infl", + "swe_out", + "eflx_lwrad_out", + "t_grnd", + "diag_out", + "qflx_evap_soi", + "eflx_soil_grnd", + "eflx_sh_tot", + "qflx_evap_veg", + "qflx_top_soil", +] + +for directory in directories: + mkdir(dir_name + "/" + directory) + +cp("$PF_SRC/test/tcl/clm/drv_clmin.dat", dir_name) +cp("$PF_SRC/test/tcl/clm/drv_vegm.dat", dir_name) +cp("$PF_SRC/test/tcl/clm/drv_vegp.dat", dir_name) +cp("$PF_SRC/test/tcl/clm/snow_forcing.1hr.txt", dir_name) + +# ----------------------------------------------------------------------------- +# File input version number +# ----------------------------------------------------------------------------- + +clm.FileVersion = 4 + +# ----------------------------------------------------------------------------- +# Process Topology +# ----------------------------------------------------------------------------- + +parser = argparse.ArgumentParser() +parser.add_argument("-p", "--p", default=1) +parser.add_argument("-q", "--q", default=1) +parser.add_argument("-r", "--r", default=1) +args = parser.parse_args() + +clm.Process.Topology.P = args.p +clm.Process.Topology.Q = args.q +clm.Process.Topology.R = args.r + +# ----------------------------------------------------------------------------- +# Computational Grid +# ----------------------------------------------------------------------------- + +clm.ComputationalGrid.Lower.X = 0.0 +clm.ComputationalGrid.Lower.Y = 0.0 +clm.ComputationalGrid.Lower.Z = 0.0 + +clm.ComputationalGrid.DX = 1000.0 +clm.ComputationalGrid.DY = 1000.0 +clm.ComputationalGrid.DZ = 0.5 + +clm.ComputationalGrid.NX = 5 +clm.ComputationalGrid.NY = 5 +clm.ComputationalGrid.NZ = 10 + +# ----------------------------------------------------------------------------- +# The Names of the GeomInputs +# ----------------------------------------------------------------------------- + +clm.GeomInput.Names = "domain_input" + +# ----------------------------------------------------------------------------- +# Domain Geometry Input +# ----------------------------------------------------------------------------- + +clm.GeomInput.domain_input.InputType = "Box" +clm.GeomInput.domain_input.GeomName = "domain" + +# ----------------------------------------------------------------------------- +# Domain Geometry +# ----------------------------------------------------------------------------- + +clm.Geom.domain.Lower.X = 0.0 +clm.Geom.domain.Lower.Y = 0.0 +clm.Geom.domain.Lower.Z = 0.0 + +clm.Geom.domain.Upper.X = 5000.0 +clm.Geom.domain.Upper.Y = 5000.0 +clm.Geom.domain.Upper.Z = 5.0 + +clm.Geom.domain.Patches = "x_lower x_upper y_lower y_upper z_lower z_upper" + +# ----------------------------------------------------------------------------- +# Perm +# ----------------------------------------------------------------------------- + +clm.Geom.Perm.Names = "domain" +clm.Geom.domain.Perm.Type = "Constant" +clm.Geom.domain.Perm.Value = 0.2 + +clm.Perm.TensorType = "TensorByGeom" +clm.Geom.Perm.TensorByGeom.Names = "domain" + +clm.Geom.domain.Perm.TensorValX = 1.0 +clm.Geom.domain.Perm.TensorValY = 1.0 +clm.Geom.domain.Perm.TensorValZ = 1.0 + +# ----------------------------------------------------------------------------- +# Specific Storage +# ----------------------------------------------------------------------------- + +clm.SpecificStorage.Type = "Constant" +clm.SpecificStorage.GeomNames = "domain" +clm.Geom.domain.SpecificStorage.Value = 1.0e-6 + +# ----------------------------------------------------------------------------- +# Phases +# ----------------------------------------------------------------------------- + +clm.Phase.Names = "water" +clm.Phase.water.Density.Type = "Constant" +clm.Phase.water.Density.Value = 1.0 +clm.Phase.water.Viscosity.Type = "Constant" +clm.Phase.water.Viscosity.Value = 1.0 + +# ----------------------------------------------------------------------------- +# Contaminants +# ----------------------------------------------------------------------------- + +clm.Contaminants.Names = "" + +# ----------------------------------------------------------------------------- +# Gravity +# ----------------------------------------------------------------------------- + +clm.Gravity = 1.0 + +# ----------------------------------------------------------------------------- +# Setup timing info +# ----------------------------------------------------------------------------- + +clm.TimingInfo.BaseUnit = 1.0 +clm.TimingInfo.StartCount = 0 +clm.TimingInfo.StartTime = 0.0 +clm.TimingInfo.StopTime = 5 +clm.TimingInfo.DumpInterval = -1 +clm.TimeStep.Type = "Constant" +clm.TimeStep.Value = 1.0 + +# ----------------------------------------------------------------------------- +# Porosity +# ----------------------------------------------------------------------------- + +clm.Geom.Porosity.GeomNames = "domain" +clm.Geom.domain.Porosity.Type = "Constant" +clm.Geom.domain.Porosity.Value = 0.390 + +# ----------------------------------------------------------------------------- +# Domain +# ----------------------------------------------------------------------------- + +clm.Domain.GeomName = "domain" + +# ----------------------------------------------------------------------------- +# Mobility +# ----------------------------------------------------------------------------- + +clm.Phase.water.Mobility.Type = "Constant" +clm.Phase.water.Mobility.Value = 1.0 + +# ----------------------------------------------------------------------------- +# Relative Permeability +# ----------------------------------------------------------------------------- + +clm.Phase.RelPerm.Type = "VanGenuchten" +clm.Phase.RelPerm.GeomNames = "domain" +clm.Geom.domain.RelPerm.Alpha = 3.5 +clm.Geom.domain.RelPerm.N = 2.0 + +# ----------------------------------------------------------------------------- +# Saturation +# ----------------------------------------------------------------------------- + +clm.Phase.Saturation.Type = "VanGenuchten" +clm.Phase.Saturation.GeomNames = "domain" +clm.Geom.domain.Saturation.Alpha = 3.5 +clm.Geom.domain.Saturation.N = 2.0 +clm.Geom.domain.Saturation.SRes = 0.01 +clm.Geom.domain.Saturation.SSat = 1.0 + +# ----------------------------------------------------------------------------- +# Wells +# ----------------------------------------------------------------------------- + +clm.Wells.Names = "" + +# ----------------------------------------------------------------------------- +# Time Cycles +# ----------------------------------------------------------------------------- + +clm.Cycle.Names = "constant" +clm.Cycle.constant.Names = "alltime" +clm.Cycle.constant.alltime.Length = 1 +clm.Cycle.constant.Repeat = -1 + +# ----------------------------------------------------------------------------- +# Boundary Conditions: Pressure +# ----------------------------------------------------------------------------- + +clm.BCPressure.PatchNames = "x_lower x_upper y_lower y_upper z_lower z_upper" + +clm.Patch.x_lower.BCPressure.Type = "FluxConst" +clm.Patch.x_lower.BCPressure.Cycle = "constant" +clm.Patch.x_lower.BCPressure.alltime.Value = 0.0 + +clm.Patch.y_lower.BCPressure.Type = "FluxConst" +clm.Patch.y_lower.BCPressure.Cycle = "constant" +clm.Patch.y_lower.BCPressure.alltime.Value = 0.0 + +clm.Patch.z_lower.BCPressure.Type = "FluxConst" +clm.Patch.z_lower.BCPressure.Cycle = "constant" +clm.Patch.z_lower.BCPressure.alltime.Value = 0.0 + +clm.Patch.x_upper.BCPressure.Type = "FluxConst" +clm.Patch.x_upper.BCPressure.Cycle = "constant" +clm.Patch.x_upper.BCPressure.alltime.Value = 0.0 + +clm.Patch.y_upper.BCPressure.Type = "FluxConst" +clm.Patch.y_upper.BCPressure.Cycle = "constant" +clm.Patch.y_upper.BCPressure.alltime.Value = 0.0 + +clm.Patch.z_upper.BCPressure.Type = "OverlandFlow" +clm.Patch.z_upper.BCPressure.Cycle = "constant" +clm.Patch.z_upper.BCPressure.alltime.Value = 0.0 + +# ----------------------------------------------------------------------------- +# Topo slopes in x-direction +# ----------------------------------------------------------------------------- + +clm.TopoSlopesX.Type = "Constant" +clm.TopoSlopesX.GeomNames = "domain" +clm.TopoSlopesX.Geom.domain.Value = -0.001 + +# ----------------------------------------------------------------------------- +# Topo slopes in y-direction +# ----------------------------------------------------------------------------- + +clm.TopoSlopesY.Type = "Constant" +clm.TopoSlopesY.GeomNames = "domain" +clm.TopoSlopesY.Geom.domain.Value = 0.001 + +# ----------------------------------------------------------------------------- +# Mannings coefficient +# ----------------------------------------------------------------------------- + +clm.Mannings.Type = "Constant" +clm.Mannings.GeomNames = "domain" +clm.Mannings.Geom.domain.Value = 5.52e-6 + +# ----------------------------------------------------------------------------- +# Phase sources +# ----------------------------------------------------------------------------- + +clm.PhaseSources.water.Type = "Constant" +clm.PhaseSources.water.GeomNames = "domain" +clm.PhaseSources.water.Geom.domain.Value = 0.0 + +# ----------------------------------------------------------------------------- +# Exact solution specification for error calculations +# ----------------------------------------------------------------------------- + +clm.KnownSolution = "NoKnownSolution" + +# ----------------------------------------------------------------------------- +# Set solver parameters +# ----------------------------------------------------------------------------- + +clm.Solver = "Richards" +clm.Solver.MaxIter = 500 + +clm.Solver.Nonlinear.MaxIter = 15 +clm.Solver.Nonlinear.ResidualTol = 1e-9 +clm.Solver.Nonlinear.EtaChoice = "EtaConstant" +clm.Solver.Nonlinear.EtaValue = 0.01 +clm.Solver.Nonlinear.UseJacobian = True +clm.Solver.Nonlinear.StepTol = 1e-20 +clm.Solver.Nonlinear.Globalization = "LineSearch" +clm.Solver.Linear.KrylovDimension = 15 +clm.Solver.Linear.MaxRestart = 2 + +clm.Solver.Linear.Preconditioner = "PFMG" +clm.Solver.PrintSubsurf = False +clm.Solver.Drop = 1e-20 +clm.Solver.AbsTol = 1e-9 + +# ----------------------------------------------------------------------------- +# CLM Settings +# ----------------------------------------------------------------------------- + +clm.Solver.LSM = "CLM" +clm.Solver.CLM.MetForcing = "1D" +clm.Solver.CLM.MetFileName = "snow_forcing.1hr.txt" +clm.Solver.CLM.MetFilePath = "." + +clm.Solver.WriteSiloCLM = False +clm.Solver.WriteSiloEvapTrans = False +clm.Solver.WriteSiloOverlandBCFlux = False +clm.Solver.PrintCLM = True + +clm.Solver.CLM.Print1dOut = False +clm.Solver.BinaryOutDir = False +clm.Solver.WriteCLMBinary = False +clm.Solver.CLM.CLMDumpInterval = 1 +clm.Solver.CLM.WriteLogs = False +clm.Solver.CLM.WriteLastRST = True +clm.Solver.CLM.DailyRST = True +clm.Solver.CLM.SingleFile = True + +# ----------------------------------------------------------------------------- +# Snow Parameterization Keys - ALL AT DEFAULTS +# These are explicitly set here to document what we're testing +# ----------------------------------------------------------------------------- + +clm.Solver.CLM.SnowPartition = "CLM" +clm.Solver.CLM.WetbulbThreshold = 274.15 +clm.Solver.CLM.ThinSnowDamping = 1.0 +clm.Solver.CLM.ThinSnowThreshold = 50.0 +clm.Solver.CLM.AlbedoScheme = "CLM" +clm.Solver.CLM.AlbedoVisNew = 0.95 +clm.Solver.CLM.AlbedoNirNew = 0.65 +clm.Solver.CLM.AlbedoMin = 0.4 +clm.Solver.CLM.AlbedoDecayVis = 0.5 +clm.Solver.CLM.AlbedoDecayNir = 0.2 +clm.Solver.CLM.AlbedoAccumA = 0.94 +clm.Solver.CLM.AlbedoThawA = 0.82 + +# ----------------------------------------------------------------------------- +# Initial conditions: water pressure +# ----------------------------------------------------------------------------- + +clm.ICPressure.Type = "HydroStaticPatch" +clm.ICPressure.GeomNames = "domain" +clm.Geom.domain.ICPressure.Value = -2.0 +clm.Geom.domain.ICPressure.RefGeom = "domain" +clm.Geom.domain.ICPressure.RefPatch = "z_upper" + +# ----------------------------------------------------------------------------- +# Run and Unload the ParFlow output files +# ----------------------------------------------------------------------------- + +correct_output_dir_name = get_absolute_path("../../correct_output/" + run_name) + +clm.run(working_directory=dir_name) + +# ----------------------------------------------------------------------------- +# Tests - Compare pressure, saturation, and CLM output +# ----------------------------------------------------------------------------- + +passed = True + +for i in range(6): + timestep = str(i).rjust(5, "0") + filename = f"/{run_name}.out.press.{timestep}.pfb" + if not pf_test_file( + dir_name + filename, + correct_output_dir_name + filename, + f"Max difference in Pressure for timestep {timestep}", + ): + passed = False + filename = f"/{run_name}.out.satur.{timestep}.pfb" + if not pf_test_file( + dir_name + filename, + correct_output_dir_name + filename, + f"Max difference in Saturation for timestep {timestep}", + ): + passed = False + + if i > 0: + filename = f"/{run_name}.out.clm_output.{timestep}.C.pfb" + if not pf_test_file( + dir_name + filename, + correct_output_dir_name + filename, + f"Max difference in CLM output for timestep {timestep}", + ): + passed = False + +rm(dir_name) + +if passed: + print(f"{run_name} : PASSED") +else: + print(f"{run_name} : FAILED") + sys.exit(1) diff --git a/test/python/clm/clm_snow_jennings.py b/test/python/clm/clm_snow_jennings.py new file mode 100644 index 000000000..ff7f1c09c --- /dev/null +++ b/test/python/clm/clm_snow_jennings.py @@ -0,0 +1,421 @@ +# ----------------------------------------------------------------------------- +# CLM Snow Jennings Test +# +# This test verifies that the Jennings (2018) rain-snow partitioning method works. +# It uses bivariate logistic regression with temperature and relative humidity: +# psnow = 1 / (1 + exp(a + b*T + g*RH)) +# +# The Jennings method accounts for both temperature and humidity effects on +# precipitation phase, which is important for accurate snow accumulation +# prediction across diverse climates. +# +# Reference: Jennings et al. (2018) Nature Communications, doi:10.1038/s41467-018-03629-7 +# +# New keys tested: +# Solver.CLM.SnowPartition = "Jennings" +# Solver.CLM.JenningsCoeffA = -10.04 +# Solver.CLM.JenningsCoeffB = 1.41 +# Solver.CLM.JenningsCoeffG = 0.09 +# ----------------------------------------------------------------------------- + +import sys +import argparse + +from parflow import Run +from parflow.tools.fs import mkdir, cp, get_absolute_path, rm +from parflow.tools.compare import pf_test_file + +run_name = "clm_snow_jennings" +clm = Run(run_name, __file__) + +# ----------------------------------------------------------------------------- +# Making output directories and copying input files +# ----------------------------------------------------------------------------- + +dir_name = get_absolute_path("test_output/" + run_name) +mkdir(dir_name) + +directories = [ + "qflx_evap_grnd", + "eflx_lh_tot", + "qflx_evap_tot", + "qflx_tran_veg", + "correct_output", + "qflx_infl", + "swe_out", + "eflx_lwrad_out", + "t_grnd", + "diag_out", + "qflx_evap_soi", + "eflx_soil_grnd", + "eflx_sh_tot", + "qflx_evap_veg", + "qflx_top_soil", +] + +for directory in directories: + mkdir(dir_name + "/" + directory) + +cp("$PF_SRC/test/tcl/clm/drv_clmin.dat", dir_name) +cp("$PF_SRC/test/tcl/clm/drv_vegm.dat", dir_name) +cp("$PF_SRC/test/tcl/clm/drv_vegp.dat", dir_name) +cp("$PF_SRC/test/tcl/clm/snow_forcing.1hr.txt", dir_name) + +# ----------------------------------------------------------------------------- +# File input version number +# ----------------------------------------------------------------------------- + +clm.FileVersion = 4 + +# ----------------------------------------------------------------------------- +# Process Topology +# ----------------------------------------------------------------------------- + +parser = argparse.ArgumentParser() +parser.add_argument("-p", "--p", default=1) +parser.add_argument("-q", "--q", default=1) +parser.add_argument("-r", "--r", default=1) +args = parser.parse_args() + +clm.Process.Topology.P = args.p +clm.Process.Topology.Q = args.q +clm.Process.Topology.R = args.r + +# ----------------------------------------------------------------------------- +# Computational Grid +# ----------------------------------------------------------------------------- + +clm.ComputationalGrid.Lower.X = 0.0 +clm.ComputationalGrid.Lower.Y = 0.0 +clm.ComputationalGrid.Lower.Z = 0.0 + +clm.ComputationalGrid.DX = 1000.0 +clm.ComputationalGrid.DY = 1000.0 +clm.ComputationalGrid.DZ = 0.5 + +clm.ComputationalGrid.NX = 5 +clm.ComputationalGrid.NY = 5 +clm.ComputationalGrid.NZ = 10 + +# ----------------------------------------------------------------------------- +# The Names of the GeomInputs +# ----------------------------------------------------------------------------- + +clm.GeomInput.Names = "domain_input" + +# ----------------------------------------------------------------------------- +# Domain Geometry Input +# ----------------------------------------------------------------------------- + +clm.GeomInput.domain_input.InputType = "Box" +clm.GeomInput.domain_input.GeomName = "domain" + +# ----------------------------------------------------------------------------- +# Domain Geometry +# ----------------------------------------------------------------------------- + +clm.Geom.domain.Lower.X = 0.0 +clm.Geom.domain.Lower.Y = 0.0 +clm.Geom.domain.Lower.Z = 0.0 + +clm.Geom.domain.Upper.X = 5000.0 +clm.Geom.domain.Upper.Y = 5000.0 +clm.Geom.domain.Upper.Z = 5.0 + +clm.Geom.domain.Patches = "x_lower x_upper y_lower y_upper z_lower z_upper" + +# ----------------------------------------------------------------------------- +# Perm +# ----------------------------------------------------------------------------- + +clm.Geom.Perm.Names = "domain" +clm.Geom.domain.Perm.Type = "Constant" +clm.Geom.domain.Perm.Value = 0.2 + +clm.Perm.TensorType = "TensorByGeom" +clm.Geom.Perm.TensorByGeom.Names = "domain" + +clm.Geom.domain.Perm.TensorValX = 1.0 +clm.Geom.domain.Perm.TensorValY = 1.0 +clm.Geom.domain.Perm.TensorValZ = 1.0 + +# ----------------------------------------------------------------------------- +# Specific Storage +# ----------------------------------------------------------------------------- + +clm.SpecificStorage.Type = "Constant" +clm.SpecificStorage.GeomNames = "domain" +clm.Geom.domain.SpecificStorage.Value = 1.0e-6 + +# ----------------------------------------------------------------------------- +# Phases +# ----------------------------------------------------------------------------- + +clm.Phase.Names = "water" +clm.Phase.water.Density.Type = "Constant" +clm.Phase.water.Density.Value = 1.0 +clm.Phase.water.Viscosity.Type = "Constant" +clm.Phase.water.Viscosity.Value = 1.0 + +# ----------------------------------------------------------------------------- +# Contaminants +# ----------------------------------------------------------------------------- + +clm.Contaminants.Names = "" + +# ----------------------------------------------------------------------------- +# Gravity +# ----------------------------------------------------------------------------- + +clm.Gravity = 1.0 + +# ----------------------------------------------------------------------------- +# Setup timing info +# ----------------------------------------------------------------------------- + +clm.TimingInfo.BaseUnit = 1.0 +clm.TimingInfo.StartCount = 0 +clm.TimingInfo.StartTime = 0.0 +clm.TimingInfo.StopTime = 5 +clm.TimingInfo.DumpInterval = -1 +clm.TimeStep.Type = "Constant" +clm.TimeStep.Value = 1.0 + +# ----------------------------------------------------------------------------- +# Porosity +# ----------------------------------------------------------------------------- + +clm.Geom.Porosity.GeomNames = "domain" +clm.Geom.domain.Porosity.Type = "Constant" +clm.Geom.domain.Porosity.Value = 0.390 + +# ----------------------------------------------------------------------------- +# Domain +# ----------------------------------------------------------------------------- + +clm.Domain.GeomName = "domain" + +# ----------------------------------------------------------------------------- +# Mobility +# ----------------------------------------------------------------------------- + +clm.Phase.water.Mobility.Type = "Constant" +clm.Phase.water.Mobility.Value = 1.0 + +# ----------------------------------------------------------------------------- +# Relative Permeability +# ----------------------------------------------------------------------------- + +clm.Phase.RelPerm.Type = "VanGenuchten" +clm.Phase.RelPerm.GeomNames = "domain" +clm.Geom.domain.RelPerm.Alpha = 3.5 +clm.Geom.domain.RelPerm.N = 2.0 + +# ----------------------------------------------------------------------------- +# Saturation +# ----------------------------------------------------------------------------- + +clm.Phase.Saturation.Type = "VanGenuchten" +clm.Phase.Saturation.GeomNames = "domain" +clm.Geom.domain.Saturation.Alpha = 3.5 +clm.Geom.domain.Saturation.N = 2.0 +clm.Geom.domain.Saturation.SRes = 0.01 +clm.Geom.domain.Saturation.SSat = 1.0 + +# ----------------------------------------------------------------------------- +# Wells +# ----------------------------------------------------------------------------- + +clm.Wells.Names = "" + +# ----------------------------------------------------------------------------- +# Time Cycles +# ----------------------------------------------------------------------------- + +clm.Cycle.Names = "constant" +clm.Cycle.constant.Names = "alltime" +clm.Cycle.constant.alltime.Length = 1 +clm.Cycle.constant.Repeat = -1 + +# ----------------------------------------------------------------------------- +# Boundary Conditions: Pressure +# ----------------------------------------------------------------------------- + +clm.BCPressure.PatchNames = "x_lower x_upper y_lower y_upper z_lower z_upper" + +clm.Patch.x_lower.BCPressure.Type = "FluxConst" +clm.Patch.x_lower.BCPressure.Cycle = "constant" +clm.Patch.x_lower.BCPressure.alltime.Value = 0.0 + +clm.Patch.y_lower.BCPressure.Type = "FluxConst" +clm.Patch.y_lower.BCPressure.Cycle = "constant" +clm.Patch.y_lower.BCPressure.alltime.Value = 0.0 + +clm.Patch.z_lower.BCPressure.Type = "FluxConst" +clm.Patch.z_lower.BCPressure.Cycle = "constant" +clm.Patch.z_lower.BCPressure.alltime.Value = 0.0 + +clm.Patch.x_upper.BCPressure.Type = "FluxConst" +clm.Patch.x_upper.BCPressure.Cycle = "constant" +clm.Patch.x_upper.BCPressure.alltime.Value = 0.0 + +clm.Patch.y_upper.BCPressure.Type = "FluxConst" +clm.Patch.y_upper.BCPressure.Cycle = "constant" +clm.Patch.y_upper.BCPressure.alltime.Value = 0.0 + +clm.Patch.z_upper.BCPressure.Type = "OverlandFlow" +clm.Patch.z_upper.BCPressure.Cycle = "constant" +clm.Patch.z_upper.BCPressure.alltime.Value = 0.0 + +# ----------------------------------------------------------------------------- +# Topo slopes in x-direction +# ----------------------------------------------------------------------------- + +clm.TopoSlopesX.Type = "Constant" +clm.TopoSlopesX.GeomNames = "domain" +clm.TopoSlopesX.Geom.domain.Value = -0.001 + +# ----------------------------------------------------------------------------- +# Topo slopes in y-direction +# ----------------------------------------------------------------------------- + +clm.TopoSlopesY.Type = "Constant" +clm.TopoSlopesY.GeomNames = "domain" +clm.TopoSlopesY.Geom.domain.Value = 0.001 + +# ----------------------------------------------------------------------------- +# Mannings coefficient +# ----------------------------------------------------------------------------- + +clm.Mannings.Type = "Constant" +clm.Mannings.GeomNames = "domain" +clm.Mannings.Geom.domain.Value = 5.52e-6 + +# ----------------------------------------------------------------------------- +# Phase sources +# ----------------------------------------------------------------------------- + +clm.PhaseSources.water.Type = "Constant" +clm.PhaseSources.water.GeomNames = "domain" +clm.PhaseSources.water.Geom.domain.Value = 0.0 + +# ----------------------------------------------------------------------------- +# Exact solution specification for error calculations +# ----------------------------------------------------------------------------- + +clm.KnownSolution = "NoKnownSolution" + +# ----------------------------------------------------------------------------- +# Set solver parameters +# ----------------------------------------------------------------------------- + +clm.Solver = "Richards" +clm.Solver.MaxIter = 500 + +clm.Solver.Nonlinear.MaxIter = 15 +clm.Solver.Nonlinear.ResidualTol = 1e-9 +clm.Solver.Nonlinear.EtaChoice = "EtaConstant" +clm.Solver.Nonlinear.EtaValue = 0.01 +clm.Solver.Nonlinear.UseJacobian = True +clm.Solver.Nonlinear.StepTol = 1e-20 +clm.Solver.Nonlinear.Globalization = "LineSearch" +clm.Solver.Linear.KrylovDimension = 15 +clm.Solver.Linear.MaxRestart = 2 + +clm.Solver.Linear.Preconditioner = "PFMG" +clm.Solver.PrintSubsurf = False +clm.Solver.Drop = 1e-20 +clm.Solver.AbsTol = 1e-9 + +# ----------------------------------------------------------------------------- +# CLM Settings +# ----------------------------------------------------------------------------- + +clm.Solver.LSM = "CLM" +clm.Solver.CLM.MetForcing = "1D" +clm.Solver.CLM.MetFileName = "snow_forcing.1hr.txt" +clm.Solver.CLM.MetFilePath = "." + +clm.Solver.WriteSiloCLM = False +clm.Solver.WriteSiloEvapTrans = False +clm.Solver.WriteSiloOverlandBCFlux = False +clm.Solver.PrintCLM = True + +clm.Solver.CLM.Print1dOut = False +clm.Solver.BinaryOutDir = False +clm.Solver.WriteCLMBinary = False +clm.Solver.CLM.CLMDumpInterval = 1 +clm.Solver.CLM.WriteLogs = False +clm.Solver.CLM.WriteLastRST = True +clm.Solver.CLM.DailyRST = True +clm.Solver.CLM.SingleFile = True + +# ----------------------------------------------------------------------------- +# Snow Partitioning - Use Jennings (2018) bivariate logistic method +# psnow = 1 / (1 + exp(a + b*T + g*RH)) +# T in Celsius, RH in percent +# These are the default coefficients from Jennings et al. (2018) +# ----------------------------------------------------------------------------- + +clm.Solver.CLM.SnowPartition = "Jennings" +clm.Solver.CLM.JenningsCoeffA = -10.04 +clm.Solver.CLM.JenningsCoeffB = 1.41 +clm.Solver.CLM.JenningsCoeffG = 0.09 + +# ----------------------------------------------------------------------------- +# Initial conditions: water pressure +# ----------------------------------------------------------------------------- + +clm.ICPressure.Type = "HydroStaticPatch" +clm.ICPressure.GeomNames = "domain" +clm.Geom.domain.ICPressure.Value = -2.0 +clm.Geom.domain.ICPressure.RefGeom = "domain" +clm.Geom.domain.ICPressure.RefPatch = "z_upper" + +# ----------------------------------------------------------------------------- +# Run ParFlow +# ----------------------------------------------------------------------------- + +clm.run(working_directory=dir_name) + +# ----------------------------------------------------------------------------- +# Tests - Compare pressure, saturation, and CLM output +# ----------------------------------------------------------------------------- + +correct_output_dir_name = get_absolute_path("$PF_SRC/test/correct_output/" + run_name) + +passed = True + +for i in range(6): + timestep = str(i).rjust(5, "0") + filename = f"/{run_name}.out.press.{timestep}.pfb" + if not pf_test_file( + dir_name + filename, + correct_output_dir_name + filename, + f"Max difference in Pressure for timestep {timestep}", + ): + passed = False + filename = f"/{run_name}.out.satur.{timestep}.pfb" + if not pf_test_file( + dir_name + filename, + correct_output_dir_name + filename, + f"Max difference in Saturation for timestep {timestep}", + ): + passed = False + + if i > 0: + filename = f"/{run_name}.out.clm_output.{timestep}.C.pfb" + if not pf_test_file( + dir_name + filename, + correct_output_dir_name + filename, + f"Max difference in CLM output for timestep {timestep}", + ): + passed = False + +rm(dir_name) + +if passed: + print(f"{run_name} : PASSED") +else: + print(f"{run_name} : FAILED") + sys.exit(1) diff --git a/test/python/clm/clm_snow_partition.py b/test/python/clm/clm_snow_partition.py new file mode 100644 index 000000000..7c809c7b9 --- /dev/null +++ b/test/python/clm/clm_snow_partition.py @@ -0,0 +1,410 @@ +# ----------------------------------------------------------------------------- +# CLM Snow Partition Test +# +# This test verifies that the rain-snow partitioning options work correctly. +# It runs with WetbulbThreshold partitioning and compares to reference output. +# +# The wet-bulb partitioning method accounts for evaporative cooling of +# falling hydrometeors, which is important in dry mountain climates where +# snow can persist at air temperatures above 0°C. +# +# Reference: Wang et al. (2019) GRL, doi:10.1029/2019GL085722 +# ----------------------------------------------------------------------------- + +import sys +import argparse + +from parflow import Run +from parflow.tools.fs import mkdir, cp, get_absolute_path, rm +from parflow.tools.compare import pf_test_file + +run_name = "clm_snow_partition" +clm = Run(run_name, __file__) + +# ----------------------------------------------------------------------------- +# Making output directories and copying input files +# ----------------------------------------------------------------------------- + +dir_name = get_absolute_path("test_output/" + run_name) +mkdir(dir_name) + +directories = [ + "qflx_evap_grnd", + "eflx_lh_tot", + "qflx_evap_tot", + "qflx_tran_veg", + "correct_output", + "qflx_infl", + "swe_out", + "eflx_lwrad_out", + "t_grnd", + "diag_out", + "qflx_evap_soi", + "eflx_soil_grnd", + "eflx_sh_tot", + "qflx_evap_veg", + "qflx_top_soil", +] + +for directory in directories: + mkdir(dir_name + "/" + directory) + +cp("$PF_SRC/test/tcl/clm/drv_clmin.dat", dir_name) +cp("$PF_SRC/test/tcl/clm/drv_vegm.dat", dir_name) +cp("$PF_SRC/test/tcl/clm/drv_vegp.dat", dir_name) +cp("$PF_SRC/test/tcl/clm/snow_forcing.1hr.txt", dir_name) + +# ----------------------------------------------------------------------------- +# File input version number +# ----------------------------------------------------------------------------- + +clm.FileVersion = 4 + +# ----------------------------------------------------------------------------- +# Process Topology +# ----------------------------------------------------------------------------- + +parser = argparse.ArgumentParser() +parser.add_argument("-p", "--p", default=1) +parser.add_argument("-q", "--q", default=1) +parser.add_argument("-r", "--r", default=1) +args = parser.parse_args() + +clm.Process.Topology.P = args.p +clm.Process.Topology.Q = args.q +clm.Process.Topology.R = args.r + +# ----------------------------------------------------------------------------- +# Computational Grid +# ----------------------------------------------------------------------------- + +clm.ComputationalGrid.Lower.X = 0.0 +clm.ComputationalGrid.Lower.Y = 0.0 +clm.ComputationalGrid.Lower.Z = 0.0 + +clm.ComputationalGrid.DX = 1000.0 +clm.ComputationalGrid.DY = 1000.0 +clm.ComputationalGrid.DZ = 0.5 + +clm.ComputationalGrid.NX = 5 +clm.ComputationalGrid.NY = 5 +clm.ComputationalGrid.NZ = 10 + +# ----------------------------------------------------------------------------- +# The Names of the GeomInputs +# ----------------------------------------------------------------------------- + +clm.GeomInput.Names = "domain_input" + +# ----------------------------------------------------------------------------- +# Domain Geometry Input +# ----------------------------------------------------------------------------- + +clm.GeomInput.domain_input.InputType = "Box" +clm.GeomInput.domain_input.GeomName = "domain" + +# ----------------------------------------------------------------------------- +# Domain Geometry +# ----------------------------------------------------------------------------- + +clm.Geom.domain.Lower.X = 0.0 +clm.Geom.domain.Lower.Y = 0.0 +clm.Geom.domain.Lower.Z = 0.0 + +clm.Geom.domain.Upper.X = 5000.0 +clm.Geom.domain.Upper.Y = 5000.0 +clm.Geom.domain.Upper.Z = 5.0 + +clm.Geom.domain.Patches = "x_lower x_upper y_lower y_upper z_lower z_upper" + +# ----------------------------------------------------------------------------- +# Perm +# ----------------------------------------------------------------------------- + +clm.Geom.Perm.Names = "domain" +clm.Geom.domain.Perm.Type = "Constant" +clm.Geom.domain.Perm.Value = 0.2 + +clm.Perm.TensorType = "TensorByGeom" +clm.Geom.Perm.TensorByGeom.Names = "domain" + +clm.Geom.domain.Perm.TensorValX = 1.0 +clm.Geom.domain.Perm.TensorValY = 1.0 +clm.Geom.domain.Perm.TensorValZ = 1.0 + +# ----------------------------------------------------------------------------- +# Specific Storage +# ----------------------------------------------------------------------------- + +clm.SpecificStorage.Type = "Constant" +clm.SpecificStorage.GeomNames = "domain" +clm.Geom.domain.SpecificStorage.Value = 1.0e-6 + +# ----------------------------------------------------------------------------- +# Phases +# ----------------------------------------------------------------------------- + +clm.Phase.Names = "water" +clm.Phase.water.Density.Type = "Constant" +clm.Phase.water.Density.Value = 1.0 +clm.Phase.water.Viscosity.Type = "Constant" +clm.Phase.water.Viscosity.Value = 1.0 + +# ----------------------------------------------------------------------------- +# Contaminants +# ----------------------------------------------------------------------------- + +clm.Contaminants.Names = "" + +# ----------------------------------------------------------------------------- +# Gravity +# ----------------------------------------------------------------------------- + +clm.Gravity = 1.0 + +# ----------------------------------------------------------------------------- +# Setup timing info +# ----------------------------------------------------------------------------- + +clm.TimingInfo.BaseUnit = 1.0 +clm.TimingInfo.StartCount = 0 +clm.TimingInfo.StartTime = 0.0 +clm.TimingInfo.StopTime = 5 +clm.TimingInfo.DumpInterval = -1 +clm.TimeStep.Type = "Constant" +clm.TimeStep.Value = 1.0 + +# ----------------------------------------------------------------------------- +# Porosity +# ----------------------------------------------------------------------------- + +clm.Geom.Porosity.GeomNames = "domain" +clm.Geom.domain.Porosity.Type = "Constant" +clm.Geom.domain.Porosity.Value = 0.390 + +# ----------------------------------------------------------------------------- +# Domain +# ----------------------------------------------------------------------------- + +clm.Domain.GeomName = "domain" + +# ----------------------------------------------------------------------------- +# Mobility +# ----------------------------------------------------------------------------- + +clm.Phase.water.Mobility.Type = "Constant" +clm.Phase.water.Mobility.Value = 1.0 + +# ----------------------------------------------------------------------------- +# Relative Permeability +# ----------------------------------------------------------------------------- + +clm.Phase.RelPerm.Type = "VanGenuchten" +clm.Phase.RelPerm.GeomNames = "domain" +clm.Geom.domain.RelPerm.Alpha = 3.5 +clm.Geom.domain.RelPerm.N = 2.0 + +# ----------------------------------------------------------------------------- +# Saturation +# ----------------------------------------------------------------------------- + +clm.Phase.Saturation.Type = "VanGenuchten" +clm.Phase.Saturation.GeomNames = "domain" +clm.Geom.domain.Saturation.Alpha = 3.5 +clm.Geom.domain.Saturation.N = 2.0 +clm.Geom.domain.Saturation.SRes = 0.01 +clm.Geom.domain.Saturation.SSat = 1.0 + +# ----------------------------------------------------------------------------- +# Wells +# ----------------------------------------------------------------------------- + +clm.Wells.Names = "" + +# ----------------------------------------------------------------------------- +# Time Cycles +# ----------------------------------------------------------------------------- + +clm.Cycle.Names = "constant" +clm.Cycle.constant.Names = "alltime" +clm.Cycle.constant.alltime.Length = 1 +clm.Cycle.constant.Repeat = -1 + +# ----------------------------------------------------------------------------- +# Boundary Conditions: Pressure +# ----------------------------------------------------------------------------- + +clm.BCPressure.PatchNames = "x_lower x_upper y_lower y_upper z_lower z_upper" + +clm.Patch.x_lower.BCPressure.Type = "FluxConst" +clm.Patch.x_lower.BCPressure.Cycle = "constant" +clm.Patch.x_lower.BCPressure.alltime.Value = 0.0 + +clm.Patch.y_lower.BCPressure.Type = "FluxConst" +clm.Patch.y_lower.BCPressure.Cycle = "constant" +clm.Patch.y_lower.BCPressure.alltime.Value = 0.0 + +clm.Patch.z_lower.BCPressure.Type = "FluxConst" +clm.Patch.z_lower.BCPressure.Cycle = "constant" +clm.Patch.z_lower.BCPressure.alltime.Value = 0.0 + +clm.Patch.x_upper.BCPressure.Type = "FluxConst" +clm.Patch.x_upper.BCPressure.Cycle = "constant" +clm.Patch.x_upper.BCPressure.alltime.Value = 0.0 + +clm.Patch.y_upper.BCPressure.Type = "FluxConst" +clm.Patch.y_upper.BCPressure.Cycle = "constant" +clm.Patch.y_upper.BCPressure.alltime.Value = 0.0 + +clm.Patch.z_upper.BCPressure.Type = "OverlandFlow" +clm.Patch.z_upper.BCPressure.Cycle = "constant" +clm.Patch.z_upper.BCPressure.alltime.Value = 0.0 + +# ----------------------------------------------------------------------------- +# Topo slopes in x-direction +# ----------------------------------------------------------------------------- + +clm.TopoSlopesX.Type = "Constant" +clm.TopoSlopesX.GeomNames = "domain" +clm.TopoSlopesX.Geom.domain.Value = -0.001 + +# ----------------------------------------------------------------------------- +# Topo slopes in y-direction +# ----------------------------------------------------------------------------- + +clm.TopoSlopesY.Type = "Constant" +clm.TopoSlopesY.GeomNames = "domain" +clm.TopoSlopesY.Geom.domain.Value = 0.001 + +# ----------------------------------------------------------------------------- +# Mannings coefficient +# ----------------------------------------------------------------------------- + +clm.Mannings.Type = "Constant" +clm.Mannings.GeomNames = "domain" +clm.Mannings.Geom.domain.Value = 5.52e-6 + +# ----------------------------------------------------------------------------- +# Phase sources +# ----------------------------------------------------------------------------- + +clm.PhaseSources.water.Type = "Constant" +clm.PhaseSources.water.GeomNames = "domain" +clm.PhaseSources.water.Geom.domain.Value = 0.0 + +# ----------------------------------------------------------------------------- +# Exact solution specification for error calculations +# ----------------------------------------------------------------------------- + +clm.KnownSolution = "NoKnownSolution" + +# ----------------------------------------------------------------------------- +# Set solver parameters +# ----------------------------------------------------------------------------- + +clm.Solver = "Richards" +clm.Solver.MaxIter = 500 + +clm.Solver.Nonlinear.MaxIter = 15 +clm.Solver.Nonlinear.ResidualTol = 1e-9 +clm.Solver.Nonlinear.EtaChoice = "EtaConstant" +clm.Solver.Nonlinear.EtaValue = 0.01 +clm.Solver.Nonlinear.UseJacobian = True +clm.Solver.Nonlinear.StepTol = 1e-20 +clm.Solver.Nonlinear.Globalization = "LineSearch" +clm.Solver.Linear.KrylovDimension = 15 +clm.Solver.Linear.MaxRestart = 2 + +clm.Solver.Linear.Preconditioner = "PFMG" +clm.Solver.PrintSubsurf = False +clm.Solver.Drop = 1e-20 +clm.Solver.AbsTol = 1e-9 + +# ----------------------------------------------------------------------------- +# CLM Settings +# ----------------------------------------------------------------------------- + +clm.Solver.LSM = "CLM" +clm.Solver.CLM.MetForcing = "1D" +clm.Solver.CLM.MetFileName = "snow_forcing.1hr.txt" +clm.Solver.CLM.MetFilePath = "." + +clm.Solver.WriteSiloCLM = False +clm.Solver.WriteSiloEvapTrans = False +clm.Solver.WriteSiloOverlandBCFlux = False +clm.Solver.PrintCLM = True + +clm.Solver.CLM.Print1dOut = False +clm.Solver.BinaryOutDir = False +clm.Solver.WriteCLMBinary = False +clm.Solver.CLM.CLMDumpInterval = 1 +clm.Solver.CLM.WriteLogs = False +clm.Solver.CLM.WriteLastRST = True +clm.Solver.CLM.DailyRST = True +clm.Solver.CLM.SingleFile = True + +# ----------------------------------------------------------------------------- +# Snow Partitioning - Use WetbulbThreshold method +# This exercises the wet-bulb temperature partitioning code path +# ----------------------------------------------------------------------------- + +clm.Solver.CLM.SnowPartition = "WetbulbThreshold" +clm.Solver.CLM.WetbulbThreshold = 274.15 # 1 degree C + +# ----------------------------------------------------------------------------- +# Initial conditions: water pressure +# ----------------------------------------------------------------------------- + +clm.ICPressure.Type = "HydroStaticPatch" +clm.ICPressure.GeomNames = "domain" +clm.Geom.domain.ICPressure.Value = -2.0 +clm.Geom.domain.ICPressure.RefGeom = "domain" +clm.Geom.domain.ICPressure.RefPatch = "z_upper" + +# ----------------------------------------------------------------------------- +# Run and Unload the ParFlow output files +# ----------------------------------------------------------------------------- + +correct_output_dir_name = get_absolute_path("../../correct_output/" + run_name) + +clm.run(working_directory=dir_name) + +# ----------------------------------------------------------------------------- +# Tests - Compare pressure, saturation, and CLM output +# ----------------------------------------------------------------------------- + +passed = True + +for i in range(6): + timestep = str(i).rjust(5, "0") + filename = f"/{run_name}.out.press.{timestep}.pfb" + if not pf_test_file( + dir_name + filename, + correct_output_dir_name + filename, + f"Max difference in Pressure for timestep {timestep}", + ): + passed = False + filename = f"/{run_name}.out.satur.{timestep}.pfb" + if not pf_test_file( + dir_name + filename, + correct_output_dir_name + filename, + f"Max difference in Saturation for timestep {timestep}", + ): + passed = False + + if i > 0: + filename = f"/{run_name}.out.clm_output.{timestep}.C.pfb" + if not pf_test_file( + dir_name + filename, + correct_output_dir_name + filename, + f"Max difference in CLM output for timestep {timestep}", + ): + passed = False + +rm(dir_name) + +if passed: + print(f"{run_name} : PASSED") +else: + print(f"{run_name} : FAILED") + sys.exit(1) diff --git a/test/python/clm/clm_snow_sza_damping.py b/test/python/clm/clm_snow_sza_damping.py new file mode 100644 index 000000000..3b1b5daae --- /dev/null +++ b/test/python/clm/clm_snow_sza_damping.py @@ -0,0 +1,420 @@ +# ----------------------------------------------------------------------------- +# CLM Snow SZA Damping Test +# +# This test verifies that the solar zenith angle (SZA) based snow melt damping +# works correctly. At high solar zenith angles (low sun), CLM underestimates +# snow albedo, causing excessive melt. This damping compensates by reducing +# melt energy at high SZA. +# +# The damping varies linearly from 1.0 (no damping) at coszen >= reference +# to the damping factor at coszen <= minimum. +# +# Reference: Dang et al. (2019) The Cryosphere, doi:10.5194/tc-13-2325-2019 +# +# New keys tested: +# Solver.CLM.SZASnowDamping = 0.8 # 20% energy reduction at high SZA +# Solver.CLM.SZADampingCoszenRef = 0.5 # Reference coszen (60 deg) +# Solver.CLM.SZADampingCoszenMin = 0.1 # Min coszen (84 deg) +# ----------------------------------------------------------------------------- + +import sys +import argparse + +from parflow import Run +from parflow.tools.fs import mkdir, cp, get_absolute_path, rm +from parflow.tools.compare import pf_test_file + +run_name = "clm_snow_sza_damping" +clm = Run(run_name, __file__) + +# ----------------------------------------------------------------------------- +# Making output directories and copying input files +# ----------------------------------------------------------------------------- + +dir_name = get_absolute_path("test_output/" + run_name) +mkdir(dir_name) + +directories = [ + "qflx_evap_grnd", + "eflx_lh_tot", + "qflx_evap_tot", + "qflx_tran_veg", + "correct_output", + "qflx_infl", + "swe_out", + "eflx_lwrad_out", + "t_grnd", + "diag_out", + "qflx_evap_soi", + "eflx_soil_grnd", + "eflx_sh_tot", + "qflx_evap_veg", + "qflx_top_soil", +] + +for directory in directories: + mkdir(dir_name + "/" + directory) + +cp("$PF_SRC/test/tcl/clm/drv_clmin.dat", dir_name) +cp("$PF_SRC/test/tcl/clm/drv_vegm.dat", dir_name) +cp("$PF_SRC/test/tcl/clm/drv_vegp.dat", dir_name) +cp("$PF_SRC/test/tcl/clm/snow_forcing.1hr.txt", dir_name) + +# ----------------------------------------------------------------------------- +# File input version number +# ----------------------------------------------------------------------------- + +clm.FileVersion = 4 + +# ----------------------------------------------------------------------------- +# Process Topology +# ----------------------------------------------------------------------------- + +parser = argparse.ArgumentParser() +parser.add_argument("-p", "--p", default=1) +parser.add_argument("-q", "--q", default=1) +parser.add_argument("-r", "--r", default=1) +args = parser.parse_args() + +clm.Process.Topology.P = args.p +clm.Process.Topology.Q = args.q +clm.Process.Topology.R = args.r + +# ----------------------------------------------------------------------------- +# Computational Grid +# ----------------------------------------------------------------------------- + +clm.ComputationalGrid.Lower.X = 0.0 +clm.ComputationalGrid.Lower.Y = 0.0 +clm.ComputationalGrid.Lower.Z = 0.0 + +clm.ComputationalGrid.DX = 1000.0 +clm.ComputationalGrid.DY = 1000.0 +clm.ComputationalGrid.DZ = 0.5 + +clm.ComputationalGrid.NX = 5 +clm.ComputationalGrid.NY = 5 +clm.ComputationalGrid.NZ = 10 + +# ----------------------------------------------------------------------------- +# The Names of the GeomInputs +# ----------------------------------------------------------------------------- + +clm.GeomInput.Names = "domain_input" + +# ----------------------------------------------------------------------------- +# Domain Geometry Input +# ----------------------------------------------------------------------------- + +clm.GeomInput.domain_input.InputType = "Box" +clm.GeomInput.domain_input.GeomName = "domain" + +# ----------------------------------------------------------------------------- +# Domain Geometry +# ----------------------------------------------------------------------------- + +clm.Geom.domain.Lower.X = 0.0 +clm.Geom.domain.Lower.Y = 0.0 +clm.Geom.domain.Lower.Z = 0.0 + +clm.Geom.domain.Upper.X = 5000.0 +clm.Geom.domain.Upper.Y = 5000.0 +clm.Geom.domain.Upper.Z = 5.0 + +clm.Geom.domain.Patches = "x_lower x_upper y_lower y_upper z_lower z_upper" + +# ----------------------------------------------------------------------------- +# Perm +# ----------------------------------------------------------------------------- + +clm.Geom.Perm.Names = "domain" +clm.Geom.domain.Perm.Type = "Constant" +clm.Geom.domain.Perm.Value = 0.2 + +clm.Perm.TensorType = "TensorByGeom" +clm.Geom.Perm.TensorByGeom.Names = "domain" + +clm.Geom.domain.Perm.TensorValX = 1.0 +clm.Geom.domain.Perm.TensorValY = 1.0 +clm.Geom.domain.Perm.TensorValZ = 1.0 + +# ----------------------------------------------------------------------------- +# Specific Storage +# ----------------------------------------------------------------------------- + +clm.SpecificStorage.Type = "Constant" +clm.SpecificStorage.GeomNames = "domain" +clm.Geom.domain.SpecificStorage.Value = 1.0e-6 + +# ----------------------------------------------------------------------------- +# Phases +# ----------------------------------------------------------------------------- + +clm.Phase.Names = "water" +clm.Phase.water.Density.Type = "Constant" +clm.Phase.water.Density.Value = 1.0 +clm.Phase.water.Viscosity.Type = "Constant" +clm.Phase.water.Viscosity.Value = 1.0 + +# ----------------------------------------------------------------------------- +# Contaminants +# ----------------------------------------------------------------------------- + +clm.Contaminants.Names = "" + +# ----------------------------------------------------------------------------- +# Gravity +# ----------------------------------------------------------------------------- + +clm.Gravity = 1.0 + +# ----------------------------------------------------------------------------- +# Setup timing info +# ----------------------------------------------------------------------------- + +clm.TimingInfo.BaseUnit = 1.0 +clm.TimingInfo.StartCount = 0 +clm.TimingInfo.StartTime = 0.0 +clm.TimingInfo.StopTime = 5 +clm.TimingInfo.DumpInterval = -1 +clm.TimeStep.Type = "Constant" +clm.TimeStep.Value = 1.0 + +# ----------------------------------------------------------------------------- +# Porosity +# ----------------------------------------------------------------------------- + +clm.Geom.Porosity.GeomNames = "domain" +clm.Geom.domain.Porosity.Type = "Constant" +clm.Geom.domain.Porosity.Value = 0.390 + +# ----------------------------------------------------------------------------- +# Domain +# ----------------------------------------------------------------------------- + +clm.Domain.GeomName = "domain" + +# ----------------------------------------------------------------------------- +# Mobility +# ----------------------------------------------------------------------------- + +clm.Phase.water.Mobility.Type = "Constant" +clm.Phase.water.Mobility.Value = 1.0 + +# ----------------------------------------------------------------------------- +# Relative Permeability +# ----------------------------------------------------------------------------- + +clm.Phase.RelPerm.Type = "VanGenuchten" +clm.Phase.RelPerm.GeomNames = "domain" +clm.Geom.domain.RelPerm.Alpha = 3.5 +clm.Geom.domain.RelPerm.N = 2.0 + +# ----------------------------------------------------------------------------- +# Saturation +# ----------------------------------------------------------------------------- + +clm.Phase.Saturation.Type = "VanGenuchten" +clm.Phase.Saturation.GeomNames = "domain" +clm.Geom.domain.Saturation.Alpha = 3.5 +clm.Geom.domain.Saturation.N = 2.0 +clm.Geom.domain.Saturation.SRes = 0.01 +clm.Geom.domain.Saturation.SSat = 1.0 + +# ----------------------------------------------------------------------------- +# Wells +# ----------------------------------------------------------------------------- + +clm.Wells.Names = "" + +# ----------------------------------------------------------------------------- +# Time Cycles +# ----------------------------------------------------------------------------- + +clm.Cycle.Names = "constant" +clm.Cycle.constant.Names = "alltime" +clm.Cycle.constant.alltime.Length = 1 +clm.Cycle.constant.Repeat = -1 + +# ----------------------------------------------------------------------------- +# Boundary Conditions: Pressure +# ----------------------------------------------------------------------------- + +clm.BCPressure.PatchNames = "x_lower x_upper y_lower y_upper z_lower z_upper" + +clm.Patch.x_lower.BCPressure.Type = "FluxConst" +clm.Patch.x_lower.BCPressure.Cycle = "constant" +clm.Patch.x_lower.BCPressure.alltime.Value = 0.0 + +clm.Patch.y_lower.BCPressure.Type = "FluxConst" +clm.Patch.y_lower.BCPressure.Cycle = "constant" +clm.Patch.y_lower.BCPressure.alltime.Value = 0.0 + +clm.Patch.z_lower.BCPressure.Type = "FluxConst" +clm.Patch.z_lower.BCPressure.Cycle = "constant" +clm.Patch.z_lower.BCPressure.alltime.Value = 0.0 + +clm.Patch.x_upper.BCPressure.Type = "FluxConst" +clm.Patch.x_upper.BCPressure.Cycle = "constant" +clm.Patch.x_upper.BCPressure.alltime.Value = 0.0 + +clm.Patch.y_upper.BCPressure.Type = "FluxConst" +clm.Patch.y_upper.BCPressure.Cycle = "constant" +clm.Patch.y_upper.BCPressure.alltime.Value = 0.0 + +clm.Patch.z_upper.BCPressure.Type = "OverlandFlow" +clm.Patch.z_upper.BCPressure.Cycle = "constant" +clm.Patch.z_upper.BCPressure.alltime.Value = 0.0 + +# ----------------------------------------------------------------------------- +# Topo slopes in x-direction +# ----------------------------------------------------------------------------- + +clm.TopoSlopesX.Type = "Constant" +clm.TopoSlopesX.GeomNames = "domain" +clm.TopoSlopesX.Geom.domain.Value = -0.001 + +# ----------------------------------------------------------------------------- +# Topo slopes in y-direction +# ----------------------------------------------------------------------------- + +clm.TopoSlopesY.Type = "Constant" +clm.TopoSlopesY.GeomNames = "domain" +clm.TopoSlopesY.Geom.domain.Value = 0.001 + +# ----------------------------------------------------------------------------- +# Mannings coefficient +# ----------------------------------------------------------------------------- + +clm.Mannings.Type = "Constant" +clm.Mannings.GeomNames = "domain" +clm.Mannings.Geom.domain.Value = 5.52e-6 + +# ----------------------------------------------------------------------------- +# Phase sources +# ----------------------------------------------------------------------------- + +clm.PhaseSources.water.Type = "Constant" +clm.PhaseSources.water.GeomNames = "domain" +clm.PhaseSources.water.Geom.domain.Value = 0.0 + +# ----------------------------------------------------------------------------- +# Exact solution specification for error calculations +# ----------------------------------------------------------------------------- + +clm.KnownSolution = "NoKnownSolution" + +# ----------------------------------------------------------------------------- +# Set solver parameters +# ----------------------------------------------------------------------------- + +clm.Solver = "Richards" +clm.Solver.MaxIter = 500 + +clm.Solver.Nonlinear.MaxIter = 15 +clm.Solver.Nonlinear.ResidualTol = 1e-9 +clm.Solver.Nonlinear.EtaChoice = "EtaConstant" +clm.Solver.Nonlinear.EtaValue = 0.01 +clm.Solver.Nonlinear.UseJacobian = True +clm.Solver.Nonlinear.StepTol = 1e-20 +clm.Solver.Nonlinear.Globalization = "LineSearch" +clm.Solver.Linear.KrylovDimension = 15 +clm.Solver.Linear.MaxRestart = 2 + +clm.Solver.Linear.Preconditioner = "PFMG" +clm.Solver.PrintSubsurf = False +clm.Solver.Drop = 1e-20 +clm.Solver.AbsTol = 1e-9 + +# ----------------------------------------------------------------------------- +# CLM Settings +# ----------------------------------------------------------------------------- + +clm.Solver.LSM = "CLM" +clm.Solver.CLM.MetForcing = "1D" +clm.Solver.CLM.MetFileName = "snow_forcing.1hr.txt" +clm.Solver.CLM.MetFilePath = "." + +clm.Solver.WriteSiloCLM = False +clm.Solver.WriteSiloEvapTrans = False +clm.Solver.WriteSiloOverlandBCFlux = False +clm.Solver.PrintCLM = True + +clm.Solver.CLM.Print1dOut = False +clm.Solver.BinaryOutDir = False +clm.Solver.WriteCLMBinary = False +clm.Solver.CLM.CLMDumpInterval = 1 +clm.Solver.CLM.WriteLogs = False +clm.Solver.CLM.WriteLastRST = True +clm.Solver.CLM.DailyRST = True +clm.Solver.CLM.SingleFile = True + +# ----------------------------------------------------------------------------- +# SZA-Based Snow Melt Damping +# At high solar zenith angles (low sun), CLM underestimates snow albedo. +# This damping reduces melt energy to compensate. +# coszen_ref = 0.5 corresponds to SZA = 60 degrees (CLM assumption) +# coszen_min = 0.1 corresponds to SZA = 84 degrees +# ----------------------------------------------------------------------------- + +clm.Solver.CLM.SZASnowDamping = 0.8 # 20% reduction at high SZA +clm.Solver.CLM.SZADampingCoszenRef = 0.5 # Start damping below SZA=60 deg +clm.Solver.CLM.SZADampingCoszenMin = 0.1 # Max damping at SZA=84 deg + +# ----------------------------------------------------------------------------- +# Initial conditions: water pressure +# ----------------------------------------------------------------------------- + +clm.ICPressure.Type = "HydroStaticPatch" +clm.ICPressure.GeomNames = "domain" +clm.Geom.domain.ICPressure.Value = -2.0 +clm.Geom.domain.ICPressure.RefGeom = "domain" +clm.Geom.domain.ICPressure.RefPatch = "z_upper" + +# ----------------------------------------------------------------------------- +# Run ParFlow +# ----------------------------------------------------------------------------- + +clm.run(working_directory=dir_name) + +# ----------------------------------------------------------------------------- +# Tests - Compare pressure, saturation, and CLM output +# ----------------------------------------------------------------------------- + +correct_output_dir_name = get_absolute_path("$PF_SRC/test/correct_output/" + run_name) + +passed = True + +for i in range(6): + timestep = str(i).rjust(5, "0") + filename = f"/{run_name}.out.press.{timestep}.pfb" + if not pf_test_file( + dir_name + filename, + correct_output_dir_name + filename, + f"Max difference in Pressure for timestep {timestep}", + ): + passed = False + filename = f"/{run_name}.out.satur.{timestep}.pfb" + if not pf_test_file( + dir_name + filename, + correct_output_dir_name + filename, + f"Max difference in Saturation for timestep {timestep}", + ): + passed = False + + if i > 0: + filename = f"/{run_name}.out.clm_output.{timestep}.C.pfb" + if not pf_test_file( + dir_name + filename, + correct_output_dir_name + filename, + f"Max difference in CLM output for timestep {timestep}", + ): + passed = False + +rm(dir_name) + +if passed: + print(f"{run_name} : PASSED") +else: + print(f"{run_name} : FAILED") + sys.exit(1) diff --git a/test/tcl/clm/snow_forcing.1hr.txt b/test/tcl/clm/snow_forcing.1hr.txt new file mode 100644 index 000000000..45fb67bc3 --- /dev/null +++ b/test/tcl/clm/snow_forcing.1hr.txt @@ -0,0 +1,6 @@ +50.0 280.0 0.0000005 274.0 1.5 -0.5 95000.0 0.0035 +25.0 275.0 0.0000008 273.0 2.0 -1.0 95100.0 0.0030 +0.0 270.0 0.0000010 271.5 1.8 -0.8 95200.0 0.0025 +0.0 268.0 0.0000012 270.0 1.5 -0.5 95300.0 0.0022 +0.0 265.0 0.0000008 269.0 1.2 -0.3 95400.0 0.0020 +25.0 270.0 0.0000005 270.5 1.0 -0.2 95500.0 0.0028