diff --git a/.vscode/launch.json b/.vscode/launch.json new file mode 100644 index 00000000..264377de --- /dev/null +++ b/.vscode/launch.json @@ -0,0 +1,35 @@ +{ + // Use IntelliSense to learn about possible attributes. + // Hover to view descriptions of existing attributes. + // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387 + "version": "0.2.0", + "configurations": [ + { + "type": "OctaveDebugger", + "request": "launch", + "name": "Octave: Execute selected m-file", + "program": "${file}", + "octave": "octave.bat --no-gui --interactive --silent", + "autoTerminate": true, + "sourceFolder": "${workspaceFolder}/src" + }, + { + "type": "OctaveDebugger", + "request": "launch", + "name": "Octave: Debug STEMMUS-SCOPE", + "program": "debug_Octave.m", + "octave": "octave.bat --no-gui --interactive --silent", + "autoTerminate": true, + "sourceFolder": "${workspaceFolder}/src" + }, + { + "type": "OctaveDebugger", + "request": "launch", + "name": "Octave: Interactive", + "program": "", + "octave": "octave.bat --no-gui --interactive --silent", + "autoTerminate": false, + "sourceFolder": "${workspaceFolder}" + }, + ] +} \ No newline at end of file diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index e94694e2..b53b92be 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -1,15 +1,17 @@ # Contributing guidelines -We welcome any kind of contributions to our software, from simple -comment or question to a full [pull -request](https://help.github.com/articles/about-pull-requests/). Please +This repository includes the MATLAB source code of the STEMMUS-SCOPE model. We welcome any +kind of contributions to our software, from simple comments or questions to a full +[pull request](https://help.github.com/articles/about-pull-requests/). Please read and follow our contributing guidelines. ## Contributing via GitHub -Note when we want to work with `STEMMUS_SCOPE` repository on a new computer for -the first time, we need to configure a few things following steps 1 to 4 below. +If you want to work with the `STEMMUS_SCOPE` repository for the first time, or on a new computer, +you need to configure a few things following steps 1 through 5 below. +
+ Steps 1 to 5 ### 1. Enable two-factor authentication It is strongly recommended using two-factor authentication. Here is the link of @@ -132,10 +134,24 @@ git clone git@github.com:EcoExtreML/STEMMUS_SCOPE.git Now a new GitHub folder `STEMMUS_SCOPE` is created in your `HOME` directory. > In this command, we clone the repository using `ssh` option. As we set the ssh -connection in [**Step 2**](#2-set-ssh-connection), this command here does not ask for our user name and -password. +connection in [**Step 2**](#2-set-ssh-connection), this command here does not +ask for our user name and password. ### 5. Collaborate using GitHub To know about the most common Git commands, follow the guides -[here](https://hackmd.io/B4v6KwsBRzG-akLDF8e5pg). \ No newline at end of file +[here](https://hackmd.io/B4v6KwsBRzG-akLDF8e5pg). +
+ +## Development of the MATLAB source of STEMMUS_SCOPE model + +To contribute to the STEMMUS_SCOPE model, you need access to the model source code that is stored in the repository [STEMMUS_SCOPE](https://github.com/EcoExtreML/STEMMUS_SCOPE). You also need a MATLAB license. MATLAB `2021a` is installed on +[Snellius]((https://servicedesk.surfsara.nl/wiki/display/WIKI/Snellius)) and [CRIB](https://crib.utwente.nl/), see [this instruction](https://pystemmusscope.readthedocs.io/en/latest/contributing_link.html#development-of-stemmus-scope-model). + +## Create an executable file of STEMMUS_SCOPE + +See the [exe readme](./exe/README.md). + +## Follow MATLAB style guidelines + +Please follow the style introduced in [MATLAB Guidelines 2.0, Richard Johnson](http://cnl.sogang.ac.kr/cnlab/lectures/programming/matlab/Richard_Johnson-MatlabStyle2_book.pdf). \ No newline at end of file diff --git a/Octave_instructions.md b/Octave_instructions.md new file mode 100644 index 00000000..5cdd8e28 --- /dev/null +++ b/Octave_instructions.md @@ -0,0 +1,48 @@ +# Using STEMMUS-SCOPE with GNU Octave +The downloads can be found here +https://octave.org/download + +TODO: installation of octave on linux +After installation, launch octave and install the following Octave packages: +`pkg install -forge netcdf` +`pkg install -forge statistics` + +For off-line installation, first, download the packages [netcdf](https://octave.sourceforge.io/netcdf/index.html), [io](https://octave.sourceforge.io/io/index.html) and [statistics](https://octave.sourceforge.io/statistics/index.html). Then, launch octave and run: + +`pkg install netcdf-1.0.16.tar.gz` +`pkg install io-2.6.4.tar.gz` +`pkg install statistics-1.4.3.tar.gz` + +### VS Code setup +Add Octave to path, e.g. for (64-bit) Windows add the following folders: +`C:\Program Files\GNU Octave\Octave-7.1.0\mingw64\bin` +`C:\Program Files\GNU Octave\Octave-7.1.0\usr\bin` + +Add the extensions +`Octave Debugger` by Paulo Silva https://marketplace.visualstudio.com/items?itemName=paulosilva.vsc-octave-debugger +This allows to run Octave in the VS Code debugger. +The debugger configurations are included in `/.vscode/launch.json` + +`Octave Hacking` by Andrew Janke https://marketplace.visualstudio.com/items?itemName=apjanke.octave-hacking +This adds syntax highlighting and formatting. + +### Running STEMMUS-SCOPE in Octave +It is possible to run STEMMUS-SCOPE from the command line with the following setup: +`octave.bat --no-gui --interactive --silent --eval "STEMMUS_SCOPE_exe('path_to_config_file')"` + +On a Unix system, use `octave` instead of `octave.bat`. +### Developing STEMMUS-SCOPE in Octave +Open the `run_Octave.m` file, either in VS Code or the Octave GUI. + +#### Octave GUI +Set the workspace to the `STEMMUS_SCOPE/src` folder, and open the `run_Octave.m` file. +Here you can set the config file that should be used, and then run the file. + +#### VS Code +While having the `STEMMUS_SCOPE` folder as the workspace, open the debugger and select `Octave: Debug STEMMUS-SCOPE`. +Start the debugger to run (and debug) the model. + +In the `run_Octave.m` file you can set the config file that should be used. + + + diff --git a/README.md b/README.md index 640531e6..b1193018 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,9 @@ # STEMMUS_SCOPE -`STEMMUS_SCOPE` serves as an integrated code of SCOPE and STEMMUS. +Integrated code of SCOPE and STEMMUS. -SCOPE is a radiative transfer and energy balance model, and STEMMUS model is a two-phase mass and heat transfer model. For more information about the coupling between these two models, please check [this reference](https://gmd.copernicus.org/articles/14/1379/2021/). +SCOPE is a radiative transfer and energy balance model, and STEMMUS model is a two-phase mass and heat transfer model. For more information about the coupling between these two models, please check [this reference](https://gmd.copernicus.org/articles/14/1379/2021/). Before running the model, you need to prepare input data and a configuration file. This can be done using the python package +[PyStemmusScope](https://pystemmusscope.readthedocs.io). Logo (by Zeng & Su, 2021) @@ -41,6 +42,3 @@ have a look at the [contribution guidelines](CONTRIBUTING.md). ## How to cite us [![RSD](https://img.shields.io/badge/rsd-ecoextreml-00a3e3.svg)](https://research-software-directory.org/projects/ecoextreml) - - -More information will follow soon. diff --git a/config_file_crib.txt b/config_file_crib.txt index c5ce39cb..0d9c9dde 100644 --- a/config_file_crib.txt +++ b/config_file_crib.txt @@ -1,7 +1,14 @@ +WorkDir=/home/jovyan/private/ SoilPropertyPath=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/SoilProperty/ -InputPath=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/ -OutputPath=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/output/ ForcingPath=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/Plumber2_data/ ForcingFileName=AU-DaS_2010-2017_OzFlux_Met.nc +directional=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/directional/ +fluspect_parameters=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/fluspect_parameters/ +leafangles=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/leafangles/ +radiationdata=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/radiationdata/ +soil_spectrum=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/soil_spectrum/ +input_data=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/input_data.xlsx InitialConditionPath=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/Initial_condition/ -DurationSize=NA +NumberOfTimeSteps=NA +InputPath= +OutputPath= \ No newline at end of file diff --git a/config_file_snellius.txt b/config_file_snellius.txt deleted file mode 100644 index 7e7f7fae..00000000 --- a/config_file_snellius.txt +++ /dev/null @@ -1,8 +0,0 @@ -SoilPropertyPath=/projects/0/einf2480/model_parameters/soil_property/ -InputPath=/scratch-shared/ecoextreml/stemmus_scope/input/ -OutputPath=/scratch-shared/ecoextreml/stemmus_scope/output/ -ForcingPath=/projects/0/einf2480/forcing/plumber2_data/ -ForcingFileName=FI-Hyy_1996-2014_FLUXNET2015_Met.nc -VegetationPropertyPath=/projects/0/einf2480/model_parameters/vegetation_property/ -InitialConditionPath=/projects/0/einf2480/soil_initialcondition/ -DurationSize=17520 \ No newline at end of file diff --git a/config_file_template.txt b/config_file_template.txt new file mode 100644 index 00000000..a78c0cb4 --- /dev/null +++ b/config_file_template.txt @@ -0,0 +1,14 @@ +WorkDir=/path_to_working_directory/ +SoilPropertyPath=/path_to_soil_property_data/ +ForcingPath=/path_to_forcing_data/ +ForcingFileName=AU-DaS_2010-2017_OzFlux_Met.nc +directional=/path_to_directional_data/ +fluspect_parameters=/path_to_fluspect_parameters_data/ +leafangles=/path_to_leafangles_data/ +radiationdata=/path_to_radiation_data/ +soil_spectrum=/path_to_soil_spectra_data/ +InitialConditionPath=/path_to_soil_initial_condition_data/ +input_data=/path_to_input_data.xlsx_file/ +NumberOfTimeSteps=NA +InputPath=will_be_created_under_WorkDir +OutputPath=will_be_created_under_WorkDir diff --git a/exe/STEMMUS_SCOPE b/exe/STEMMUS_SCOPE deleted file mode 100755 index cd97ea73..00000000 Binary files a/exe/STEMMUS_SCOPE and /dev/null differ diff --git a/pull_request_template.md b/pull_request_template.md index 8aa1b1a5..de5657b4 100644 --- a/pull_request_template.md +++ b/pull_request_template.md @@ -5,9 +5,9 @@ Please add a description of the changes proposed in the pull request. Please re-generate exe file if matlab codes are changed. About how to create an exe file, see [exe readme](https://github.com/EcoExtreML/STEMMUS_SCOPE/blob/main/exe/README.md). ### Checklist -- [ ] Add a refernce to related issues. +- [ ] Add a reference to related issues. - [ ] @mentions of the person or team responsible for reviewing proposed changes. - [ ] This pull request has a descriptive title. -- [ ] Code is written according to the [guidelines].(http://cnl.sogang.ac.kr/cnlab/lectures/programming/matlab/Richard_Johnson-MatlabStyle2_book.pdf) +- [ ] Code is written according to the [guidelines](http://cnl.sogang.ac.kr/cnlab/lectures/programming/matlab/Richard_Johnson-MatlabStyle2_book.pdf). - [ ] Documentation is available. - [ ] Model runs successfully. diff --git a/run_model_on_snellius/README.md b/run_model_on_snellius/README.md new file mode 100644 index 00000000..92cc0bd6 --- /dev/null +++ b/run_model_on_snellius/README.md @@ -0,0 +1,9 @@ +## Run STEMMUS_SCOPE on Snellius + +This folder contains files that are needed for running `STEMMUS_SCOPE` on Snellius: + +- `run_stemmus_scope_snellius.sh` is the batch script for job submission +- `run_model.py` is the python script for executing the workflow +- `config_file_snellius.txt` is the model config file + +Before submitting the job via `sbatch run_stemmus_scope_snellius.sh`, make sure that conda environment `pystemmusscope` is created, see the [environment file](https://github.com/EcoExtreML/STEMMUS_SCOPE_Processing/blob/main/environment.yml). Also, set the `WorkDir` and `NumberOfTimeSteps` in the model config file e.g. `config_file_snellius.txt`. diff --git a/run_model_on_snellius/config_file_snellius.txt b/run_model_on_snellius/config_file_snellius.txt new file mode 100644 index 00000000..f1bcdb3a --- /dev/null +++ b/run_model_on_snellius/config_file_snellius.txt @@ -0,0 +1,14 @@ +WorkDir=/scratch-shared/ecoextreml/stemmus_scope/ +SoilPropertyPath=/projects/0/einf2480/model_parameters/soil_property/ +ForcingPath=/projects/0/einf2480/forcing/plumber2_data/ +ForcingFileName=FI-Hyy_1996-2014_FLUXNET2015_Met.nc +directional=/projects/0/einf2480/model_parameters/vegetation_property/directional/ +fluspect_parameters=/projects/0/einf2480/model_parameters/vegetation_property/fluspect_parameters/ +leafangles=/projects/0/einf2480/model_parameters/vegetation_property/leafangles/ +radiationdata=/projects/0/einf2480/model_parameters/vegetation_property/radiationdata/ +soil_spectrum=/projects/0/einf2480/model_parameters/vegetation_property/soil_spectrum/ +input_data=/projects/0/einf2480/model_parameters/vegetation_property/input_data.xlsx +InitialConditionPath=/projects/0/einf2480/soil_initialcondition/ +NumberOfTimeSteps=17520 +InputPath= +OutputPath= \ No newline at end of file diff --git a/exe/README.md b/run_model_on_snellius/exe/README.md similarity index 100% rename from exe/README.md rename to run_model_on_snellius/exe/README.md diff --git a/run_model_on_snellius/exe/STEMMUS_SCOPE b/run_model_on_snellius/exe/STEMMUS_SCOPE new file mode 100755 index 00000000..112ffb55 Binary files /dev/null and b/run_model_on_snellius/exe/STEMMUS_SCOPE differ diff --git a/exe/build_stemmus_scope_exe.sh b/run_model_on_snellius/exe/build_stemmus_scope_exe.sh similarity index 100% rename from exe/build_stemmus_scope_exe.sh rename to run_model_on_snellius/exe/build_stemmus_scope_exe.sh diff --git a/run_model_on_snellius/run_model.py b/run_model_on_snellius/run_model.py new file mode 100644 index 00000000..6306216c --- /dev/null +++ b/run_model_on_snellius/run_model.py @@ -0,0 +1,89 @@ +"""Run STEMMUS_SCOPE on Snellius. +This script is used in a job submission file to loop over forcing files and run +StemmusScope model on a HPC cluster using slurm workload manager. In order to +run the model, PyStemmusScope should be installed, see +https://pystemmusscope.readthedocs.io/. +""" + +import argparse +from pathlib import Path +from PyStemmusScope import StemmusScope +from PyStemmusScope import save + + +def run_model(ncfile_index, job_id): + """Workflow executer. + + Run StemmusScope model with PyStemmusScope. + """ + path_to_config_file = "./config_file_snellius.txt" + path_to_exe_file = "./exe/STEMMUS_SCOPE" + + # create an instance of the model + model = StemmusScope(config_file=path_to_config_file, exe_file=path_to_exe_file) + + # get the forcing file + forcing_filenames_list = [ + file.name for file in Path(model.config["ForcingPath"]).iterdir() + ] + # note that index starts from 1 from bash input, so ncfile_index-1 + nc_file = forcing_filenames_list[ncfile_index - 1] + station_name = nc_file.split("_")[0] + + # feed model with the correct forcing file + config_path = model.setup(ForcingFileName=nc_file) + + # run model + model_log = model.run() + + # save output in netcdf format + required_netcdf_variables = "../utils/csv_to_nc/required_netcdf_variables.csv" + netcdf_file_name = save.to_netcdf(config_path, required_netcdf_variables) + + slurm_log_msg = [ + f"Input path is {model.config['InputPath']}", + f"Output path is {model.config['OutputPath']}", + f"model log is {model_log}", + f"model outputs are in {netcdf_file_name}", + ] + + slurm_file_name = ( + Path(model.config["OutputPath"]) + / f"slurm_{job_id}_{ncfile_index}_{station_name}.out" + ) + + # create slurm log + slurm_log(slurm_file_name, slurm_log_msg) + + +def slurm_log(slurm_file_name, slurm_log_msg): + """Create slurm log file for the submitted job.""" + with open(slurm_file_name, "w+", encoding="utf-8") as file: + for line in slurm_log_msg: + file.write(line) + file.write("\n") + + +if __name__ == "__main__": + parser = argparse.ArgumentParser() + + parser.add_argument( + "-n", + "--ncfile_index", + required=True, + type=int, + default=0, + help="index of netCDF file", + ) + parser.add_argument( + "-j", + "--job_id", + required=True, + type=int, + default=0, + help="slurm job id from snellius", + ) + # get arguments + args = parser.parse_args() + + run_model(args.ncfile_index, args.job_id) diff --git a/run_model_on_snellius/run_stemmus_scope_snellius.sh b/run_model_on_snellius/run_stemmus_scope_snellius.sh new file mode 100644 index 00000000..50309996 --- /dev/null +++ b/run_model_on_snellius/run_stemmus_scope_snellius.sh @@ -0,0 +1,48 @@ +#!/bin/bash +# This is a batch script for Snellius at Surf +# usage: cd STEMMUS_SCOPE/run_model_on_snellius; sbatch run_stemmus_scope_snellius.sh + +# SLURM settings +#SBATCH -J stemmus_scope +#SBATCH -t 00:10:00 +#SBATCH --nodes=1 +#SBATCH --ntasks=32 +#SBATCH -p thin +#SBATCH --output=./slurm_%j.out +#SBATCH --error=./slurm_%j.out + +### 1. Load module needed to run the model (no need for license) +module load 2021 + +### This is for Matlab +module load MCR/R2021a.3 + +### python environment pystemmusscope is needed +### see https://pystemmusscope.readthedocs.io +. ~/mamba/bin/activate pystemmusscope + +### 2. Some security: stop script on error and undefined variables +set -euo pipefail + +### 3. Run parallel loop +# number of cores on a shared node +ncores=32 + +# number of forcing files +nfiles=170 + +# some indices for loop +k=0 +i=0 + +for k in `seq 0 5`; do + for j in `seq 1 $ncores` ; do + ( + i=$(( ncores * k + j )) + if [[ $i -le $nfiles ]]; then + python run_model.py -n $i -j ${SLURM_JOB_ID} + fi + )& + done + wait +done \ No newline at end of file diff --git a/run_stemmus_scope_snellius.sh b/run_stemmus_scope_snellius.sh deleted file mode 100644 index 84907798..00000000 --- a/run_stemmus_scope_snellius.sh +++ /dev/null @@ -1,111 +0,0 @@ -#!/bin/bash -# This is a batch script for Snellius at Surf -# usage: cd STEMMUS_SCOPE; mkdir -p slurm; sbatch run_stemmus_scope_snellius.sh - -# SLURM settings -#SBATCH -J stemmus_scope -#SBATCH -t 00:10:00 -#SBATCH --nodes=1 -#SBATCH --ntasks=32 -#SBATCH -p thin -#SBATCH --output=./slurm/slurm_%j.out -#SBATCH --error=./slurm/slurm_%j.out - -### 1. Load module needed to run the model (no need for license) -module load 2021 - -### This is for Matlab -module load MCR/R2021a.3 - -### python environment stemmus is needed to convert csv files to nc files -### see utils/csv_to_nc/README.md -. ~/mamba/bin/activate stemmus - -### 2. Some security: stop script on error and undefined variables -set -euo pipefail - -### 3. Create a function to loop over -loop_func() { - - ### 3.2 Get paths from a config file - config="config_file_snellius.txt" - source $config - - ### 3.3 loop over forcing file using i that is the function input argument - i=$1 - ncfiles=("$(echo "$ForcingPath" | tr -d '\r')"*.nc) - - ### -1 is needed because ncfiles array starts from 0 - file=${ncfiles[$i-1]} - - start_time=$(date +%s) - base_name=$(basename ${file}) - station_name="${base_name%%_*}" - timestamp=$(date +"%Y-%m-%d-%H%M") - - ### 3.4 create input directories, - work_dir="$(echo "$InputPath" | tr -d '\r')${station_name}_${timestamp}/" - mkdir -p $work_dir - - ### 3.5 copy files to work_dir, - vegetation_path="$(echo "$VegetationPropertyPath" | tr -d '\r')" - cp -r $vegetation_path/* $work_dir - - ### 3.6 update config file for ForcingFileName and InputPath, - ### !due to "/" in work_dir, the sed uses "|" instead of "/" - station_config="${work_dir}${station_name}_${timestamp}_config.txt" - sed -e "s/ForcingFileName=.*$/ForcingFileName=$base_name/g" \ - -e "s|InputPath=.*$|InputPath=$work_dir|g" <$config >$station_config - - ### 3.7 Add info to a std_out file - std_out="./slurm/slurm_${SLURM_JOB_ID}_${i}_${station_name}.out" - echo " - SoilPropertyPath is $SoilPropertyPath - InputPath is $InputPath - OutputPath is $OutputPath - ForcingPath is $ForcingPath - VegetationPropertyPath is $VegetationPropertyPath - WorkDir is $work_dir - ForcingFileName is $base_name" > $std_out - - ### 3.8 run the model and store the time and model log info - ### set matlab log dir to slurm, otherwise java.log files are created in user home dir - export MATLAB_LOG_DIR=./slurm - exe/STEMMUS_SCOPE $station_config >> $std_out - - end_time=$(date +%s) - run_time=$(expr $end_time - $start_time) - - ## 3.9 Add some information to slurm*.out file later will be used. - completed="COMPLETED" - echo "Run is $completed. Model run time is $run_time s." >> $std_out - - ## 3.10 Convert csv files to a nc file, if run is completed - if [[ -v completed ]]; - then - python utils/csv_to_nc/generate_netcdf_files.py --config_file $station_config --variable_file utils/csv_to_nc/Variables_will_be_in_NetCDF_file.csv - fi -} - -### 5. Run parallel loop -# number of cores on a shared node -ncores=32 - -# number of forcing files -nfiles=170 - -# some indices for loop -k=0 -i=0 - -for k in `seq 0 5`; do - for j in `seq 1 $ncores` ; do - ( - i=$(( ncores * k + j )) - if [[ $i -le $nfiles ]]; then - loop_func $i - fi - )& - done - wait -done \ No newline at end of file diff --git a/src/+io/bin_to_csv.m b/src/+io/bin_to_csv.m index 70b6bd0c..43575199 100644 --- a/src/+io/bin_to_csv.m +++ b/src/+io/bin_to_csv.m @@ -1,138 +1,142 @@ function bin_to_csv(fnames, V, vmax, n_col, ns, options, SoilLayer) -%% flu -flu_names = {'simulation_number','nu_iterations','year','DoY',... - 'Rntot','lEtot','Htot', ... - 'Rnctot','lEctot','Hctot','Actot', ... - 'Rnstot','lEstot','Hstot','Gtot','Resp', ... - 'aPAR','aPAR_Cab','aPAR/rad.PAR','aPAR_Wm2','PAR','rad.Eoutf','rad.Eoutf./fluxes.aPAR_Wm2','Trap','Evap','ET','GPP','NEE'}; -flu_units = {'','','','', ... - 'W m-2','W m-2','W m-2',... - 'W m-2','W m-2','W m-2','umol m-2 s-1',... - 'W m-2',' W m-2','W m-2','W m-2','umol m-2 s-1',... - 'umol m-2 s-1',' umol m-2 s-1','umol m-2 s-1','W m-2','umol m-2 s-1','W m-2','W m-2','cm s-1','cm s-1','cm s-1','Kg m-2 s-1','Kg m-2 s-1'}; -write_output(flu_names, flu_units, fnames.flu_file, n_col.flu, ns) - -%% surftemp -surftemp_names = {'simulation_number','year','DoY',... - 'Ta','Ts(1)','Ts(2)','Tcave','Tsave'}; -surftemp_units = {'','','', ... - 'C','C','C','C','C'}; -write_output(surftemp_names, surftemp_units, fnames.surftemp_file, n_col.surftemp, ns) - -%% aerodyn -aerodyn_names = {'simulation_number','year','DoY',... - 'raa','rawc','raws','ustar','rac','ras'}; -aerodyn_units = {'','','', ... - 's m-1','s m-1','s m-1', 'm s-1','s m-1','s m-1'}; -write_output(aerodyn_names, aerodyn_units, fnames.aerodyn_file, n_col.aerodyn, ns) - -%% radiation -radiation_names = {'simulation_number','year','DoY',... - 'Rin','Rli','HemisOutShort','HemisOutLong','HemisOutTot','Netshort','Netlong','Rntot'}; -radiation_units = {'','','', ... - 'W m-2','W m-2','W m-2', 'W m-2','W m-2','W m-2','W m-2','W m-2'}; -write_output(radiation_names, radiation_units, fnames.radiation_file, n_col.radiation, ns) - -%% Get soil layer information -depth = string(SoilLayer.depth); -thickness = string(SoilLayer.thickness); - -%% Prepare soil moisture data -Sim_Theta_names = [depth; thickness]; -Sim_Theta_units = repelem({'m-3 m-3'}, length(depth)); -write_output(Sim_Theta_names, Sim_Theta_units, ... -fnames.Sim_Theta_file, n_col.Sim_Theta, ns, true) - -%% Prepare soil temperature data -Sim_Temp_names = [depth; thickness]; -Sim_Temp_units = repelem({'oC'}, length(depth)); -write_output(Sim_Temp_names, Sim_Temp_units, ... -fnames.Sim_Temp_file, n_col.Sim_Temp, ns, true) - -%% spectrum (added on 19 September 2008) - -spectrum_hemis_optical_names = {'hemispherically integrated radiation spectrum'}; -spectrum_hemis_optical_units = {'W m-2 um-1'}; -write_output(spectrum_hemis_optical_names, spectrum_hemis_optical_units, fnames.spectrum_hemis_optical_file, n_col.spectrum_hemis_optical, ns, true) - -spectrum_obsdir_optical_names = {'radiance spectrum in observation direction'}; -spectrum_obsdir_optical_units = {'W m-2 sr-1 um-1'}; -write_output(spectrum_obsdir_optical_names, spectrum_obsdir_optical_units, fnames.spectrum_obsdir_optical_file, n_col.spectrum_obsdir_optical, ns, true) - -if options.calc_ebal - write_output({'thermal BlackBody emission spectrum in observation direction'}, {'W m-2 sr-1 um-1'}, ... - fnames.spectrum_obsdir_BlackBody_file, n_col.spectrum_obsdir_BlackBody, ns, true) - if options.calc_planck - write_output({'thermal emission spectrum in hemispherical direction'}, {'W m-2 sr-1 um-1'}, ... - fnames.spectrum_hemis_thermal_file, n_col.spectrum_hemis_thermal, ns, true) - write_output({'thermal emission spectrum in observation direction'}, {'W m-2 sr-1 um-1'}, ... - fnames.spectrum_obsdir_thermal_file, n_col.spectrum_obsdir_thermal, ns, true) - end -end - write_output({'irradiance'}, {'W m-2 um-1'}, ... - fnames.irradiance_spectra_file, n_col.irradiance_spectra, ns, true) - write_output({'reflectance'}, {'fraction of radiation in observation direction *pi / irradiance'}, ... - fnames.reflectance_file, n_col.reflectance, ns, true) -%% input and parameter values (added June 2012) -%write_output(fnames.pars_and_input_file, true) -%write_output(fnames.pars_and_input_short_file, true) -%% Optional Output -if options.calc_vert_profiles - write_output({'Fraction leaves in the sun, fraction of observed, fraction of observed&visible per layer'}, {'rows: simulations or time steps, columns: layer numbers'}, ... - fnames.gap_file, n_col.gap, ns, true) - write_output({'aPAR per leaf layer'}, {'umol m-2 s-1'}, ... - fnames.layer_aPAR_file, n_col.layer_aPAR, ns, true) - write_output({'aPAR by Cab per leaf layer'}, {'umol m-2 s-1'}, ... - fnames.layer_aPAR_Cab_file, n_col.layer_aPAR_Cab, ns, true) + %% flu + flu_names = {'simulation_number','nu_iterations','year','DoY',... + 'Rntot','lEtot','Htot', ... + 'Rnctot','lEctot','Hctot','Actot', ... + 'Rnstot','lEstot','Hstot','Gtot','Resp', ... + 'aPAR','aPAR_Cab','aPAR/rad.PAR','aPAR_Wm2','PAR','rad.Eoutf','rad.Eoutf./fluxes.aPAR_Wm2','Trap','Evap','ET','GPP','NEE'}; + flu_units = {'','','','', ... + 'W m-2','W m-2','W m-2',... + 'W m-2','W m-2','W m-2','umol m-2 s-1',... + 'W m-2',' W m-2','W m-2','W m-2','umol m-2 s-1',... + 'umol m-2 s-1',' umol m-2 s-1','umol m-2 s-1','W m-2','umol m-2 s-1','W m-2','W m-2','cm s-1','cm s-1','cm s-1','Kg m-2 s-1','Kg m-2 s-1'}; + write_output(flu_names, flu_units, fnames.flu_file, n_col.flu, ns) + + %% surftemp + surftemp_names = {'simulation_number','year','DoY',... + 'Ta','Ts(1)','Ts(2)','Tcave','Tsave'}; + surftemp_units = {'','','', ... + 'C','C','C','C','C'}; + write_output(surftemp_names, surftemp_units, fnames.surftemp_file, n_col.surftemp, ns) + + %% aerodyn + aerodyn_names = {'simulation_number','year','DoY',... + 'raa','rawc','raws','ustar','rac','ras'}; + aerodyn_units = {'','','', ... + 's m-1','s m-1','s m-1', 'm s-1','s m-1','s m-1'}; + write_output(aerodyn_names, aerodyn_units, fnames.aerodyn_file, n_col.aerodyn, ns) + + %% radiation + radiation_names = {'simulation_number','year','DoY',... + 'Rin','Rli','HemisOutShort','HemisOutLong','HemisOutTot','Netshort','Netlong','Rntot'}; + radiation_units = {'','','', ... + 'W m-2','W m-2','W m-2', 'W m-2','W m-2','W m-2','W m-2','W m-2'}; + write_output(radiation_names, radiation_units, fnames.radiation_file, n_col.radiation, ns) + + %% Get soil layer information + depth = arrayfun(@(x) num2str(x), SoilLayer.depth, "UniformOutput", false); + thickness = arrayfun(@(x) num2str(x), SoilLayer.thickness, "UniformOutput", false); + + %% Prepare soil moisture data + Sim_Theta_names = [depth; thickness]; + Sim_Theta_units = repelem({'m-3 m-3'}, length(depth)); + write_output(Sim_Theta_names, Sim_Theta_units, ... + fnames.Sim_Theta_file, n_col.Sim_Theta, ns, true) + + %% Prepare soil temperature data + Sim_Temp_names = [depth; thickness]; + Sim_Temp_units = repelem({'oC'}, length(depth)); + write_output(Sim_Temp_names, Sim_Temp_units, ... + fnames.Sim_Temp_file, n_col.Sim_Temp, ns, true) + + %% spectrum (added on 19 September 2008) + + spectrum_hemis_optical_names = {'hemispherically integrated radiation spectrum'}; + spectrum_hemis_optical_units = {'W m-2 um-1'}; + write_output(spectrum_hemis_optical_names, spectrum_hemis_optical_units, fnames.spectrum_hemis_optical_file, n_col.spectrum_hemis_optical, ns, true) + + spectrum_obsdir_optical_names = {'radiance spectrum in observation direction'}; + spectrum_obsdir_optical_units = {'W m-2 sr-1 um-1'}; + write_output(spectrum_obsdir_optical_names, spectrum_obsdir_optical_units, fnames.spectrum_obsdir_optical_file, n_col.spectrum_obsdir_optical, ns, true) + if options.calc_ebal - write_output({'leaf temperature of sunlit leaves, shaded leaves, and weighted average leaf temperature per layer'}, {'^oC ^oC ^oC'}, ... - fnames.leaftemp_file, n_col.leaftemp, ns, true) - write_output({'sensible heat flux per layer'}, {'W m-2'}, ... - fnames.layer_H_file, n_col.layer_H, ns, true) - write_output({'latent heat flux per layer'}, {'W m-2'}, ... - fnames.layer_LE_file, n_col.layer_LE, ns, true) - write_output({'photosynthesis per layer'}, {'umol m-2 s-1'}, ... - fnames.layer_A_file, n_col.layer_A, ns, true) - write_output({'average NPQ = 1-(fm-fo)/(fm0-fo0), per layer'}, {''}, ... - fnames.layer_NPQ_file, n_col.layer_NPQ, ns, true) - write_output({'net radiation per leaf layer'}, {'W m-2'}, ... - fnames.layer_Rn_file, n_col.layer_Rn, ns, true) + write_output({'thermal BlackBody emission spectrum in observation direction'}, {'W m-2 sr-1 um-1'}, ... + fnames.spectrum_obsdir_BlackBody_file, n_col.spectrum_obsdir_BlackBody, ns, true) + if options.calc_planck + write_output({'thermal emission spectrum in hemispherical direction'}, {'W m-2 sr-1 um-1'}, ... + fnames.spectrum_hemis_thermal_file, n_col.spectrum_hemis_thermal, ns, true) + write_output({'thermal emission spectrum in observation direction'}, {'W m-2 sr-1 um-1'}, ... + fnames.spectrum_obsdir_thermal_file, n_col.spectrum_obsdir_thermal, ns, true) + end end -if options.calc_fluor -fluorescence_names = {'supward fluorescence per layer'}; -fluorescence_units = {'W m-2'}; -write_output(fluorescence_names, fluorescence_units, fnames.layer_fluorescence_file, n_col.layer_fluorescence, ns, true) -end -end -if options.calc_fluor - write_output({'fluorescence per simulation for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1 sr-1'}, ... - fnames.fluorescence_file, n_col.fluorescence, ns, true) - if options.calc_PSI - write_output({'fluorescence per simulation for wavelengths of 640 to 850 nm, with 1 nm resolution, for PSI only'}, {'W m-2 um-1 sr-1'}, ... - fnames.fluorescencePSI_file, n_col.fluorescencePSI, ns, true) - write_output({'fluorescence per simulation for wavelengths of 640 to 850 nm, with 1 nm resolution, for PSII only'}, {'W m-2 um-1 sr-1'}, ... - fnames.fluorescencePSII_file, n_col.fluorescencePSII, ns, true) + write_output({'irradiance'}, {'W m-2 um-1'}, ... + fnames.irradiance_spectra_file, n_col.irradiance_spectra, ns, true) + write_output({'reflectance'}, {'fraction of radiation in observation direction *pi / irradiance'}, ... + fnames.reflectance_file, n_col.reflectance, ns, true) + %% input and parameter values (added June 2012) + %write_output(fnames.pars_and_input_file, true) + %write_output(fnames.pars_and_input_short_file, true) + %% Optional Output + if options.calc_vert_profiles + write_output({'Fraction leaves in the sun, fraction of observed, fraction of observed&visible per layer'}, {'rows: simulations or time steps, columns: layer numbers'}, ... + fnames.gap_file, n_col.gap, ns, true) + write_output({'aPAR per leaf layer'}, {'umol m-2 s-1'}, ... + fnames.layer_aPAR_file, n_col.layer_aPAR, ns, true) + write_output({'aPAR by Cab per leaf layer'}, {'umol m-2 s-1'}, ... + fnames.layer_aPAR_Cab_file, n_col.layer_aPAR_Cab, ns, true) + if options.calc_ebal + write_output({'leaf temperature of sunlit leaves, shaded leaves, and weighted average leaf temperature per layer'}, {'^oC ^oC ^oC'}, ... + fnames.leaftemp_file, n_col.leaftemp, ns, true) + write_output({'sensible heat flux per layer'}, {'W m-2'}, ... + fnames.layer_H_file, n_col.layer_H, ns, true) + write_output({'latent heat flux per layer'}, {'W m-2'}, ... + fnames.layer_LE_file, n_col.layer_LE, ns, true) + write_output({'photosynthesis per layer'}, {'umol m-2 s-1'}, ... + fnames.layer_A_file, n_col.layer_A, ns, true) + write_output({'average NPQ = 1-(fm-fo)/(fm0-fo0), per layer'}, {''}, ... + fnames.layer_NPQ_file, n_col.layer_NPQ, ns, true) + write_output({'net radiation per leaf layer'}, {'W m-2'}, ... + fnames.layer_Rn_file, n_col.layer_Rn, ns, true) + end + if options.calc_fluor + fluorescence_names = {'supward fluorescence per layer'}; + fluorescence_units = {'W m-2'}; + write_output(fluorescence_names, fluorescence_units, fnames.layer_fluorescence_file, n_col.layer_fluorescence, ns, true) end - write_output({'hemispherically integrated fluorescence per simulation for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1'}, ... - fnames.fluorescence_hemis_file, n_col.fluorescence_hemis, ns, true) - write_output({'total emitted fluorescence by all leaves for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1'}, ... - fnames.fluorescence_emitted_by_all_leaves_file, n_col.fluorescence_emitted_by_all_leaves, ns, true) - write_output({'total emitted fluorescence by all photosystems for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1'}, ... - fnames.fluorescence_emitted_by_all_photosystems_file, n_col.fluorescence_emitted_by_all_photosystems, ns, true) - write_output({'TOC fluorescence contribution from sunlit leaves for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1 sr-1'}, ... - fnames.fluorescence_sunlit_file, n_col.fluorescence_sunlit, ns, true) - write_output({'TOC fluorescence contribution from shaded leaves for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1 sr-1'}, ... - fnames.fluorescence_shaded_file, n_col.fluorescence_shaded, ns, true) - write_output({'TOC fluorescence contribution from from leaves and soil after scattering for wavelenghts of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1 sr-1'}, ... - fnames.fluorescence_scattered_file, n_col.fluorescence_scattered, ns, true) -end -write_output({'Bottom of canopy irradiance in the shaded fraction, and average BOC irradiance'}, {'First 2162 columns: shaded fraction. Last 2162 columns: average BOC irradiance. Unit: Wm-2 um-1'}, ... -fnames.BOC_irradiance_file, n_col.BOC_irradiance, ns, true) + end + if options.calc_fluor + write_output({'fluorescence per simulation for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1 sr-1'}, ... + fnames.fluorescence_file, n_col.fluorescence, ns, true) + if options.calc_PSI + write_output({'fluorescence per simulation for wavelengths of 640 to 850 nm, with 1 nm resolution, for PSI only'}, {'W m-2 um-1 sr-1'}, ... + fnames.fluorescencePSI_file, n_col.fluorescencePSI, ns, true) + write_output({'fluorescence per simulation for wavelengths of 640 to 850 nm, with 1 nm resolution, for PSII only'}, {'W m-2 um-1 sr-1'}, ... + fnames.fluorescencePSII_file, n_col.fluorescencePSII, ns, true) + end + write_output({'hemispherically integrated fluorescence per simulation for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1'}, ... + fnames.fluorescence_hemis_file, n_col.fluorescence_hemis, ns, true) + write_output({'total emitted fluorescence by all leaves for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1'}, ... + fnames.fluorescence_emitted_by_all_leaves_file, n_col.fluorescence_emitted_by_all_leaves, ns, true) + write_output({'total emitted fluorescence by all photosystems for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1'}, ... + fnames.fluorescence_emitted_by_all_photosystems_file, n_col.fluorescence_emitted_by_all_photosystems, ns, true) + write_output({'TOC fluorescence contribution from sunlit leaves for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1 sr-1'}, ... + fnames.fluorescence_sunlit_file, n_col.fluorescence_sunlit, ns, true) + write_output({'TOC fluorescence contribution from shaded leaves for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1 sr-1'}, ... + fnames.fluorescence_shaded_file, n_col.fluorescence_shaded, ns, true) + write_output({'TOC fluorescence contribution from from leaves and soil after scattering for wavelenghts of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1 sr-1'}, ... + fnames.fluorescence_scattered_file, n_col.fluorescence_scattered, ns, true) + end + write_output({'Bottom of canopy irradiance in the shaded fraction, and average BOC irradiance'}, {'First 2162 columns: shaded fraction. Last 2162 columns: average BOC irradiance. Unit: Wm-2 um-1'}, ... + fnames.BOC_irradiance_file, n_col.BOC_irradiance, ns, true) + + fclose('all'); -fclose('all'); + %% Delete all (temporary) .bin files + fn = fieldnames(fnames); + for k=1:numel(fn) + delete(fnames.(fn{k})) + end -%% deleting .bin -structfun(@delete, fnames) end function write_output(header, units, bin_path, f_n_col, ns, not_header) @@ -140,30 +144,35 @@ function write_output(header, units, bin_path, f_n_col, ns, not_header) not_header = false; end n_csv = strrep(bin_path, '.bin', '.csv'); - + f_csv = fopen(n_csv, 'w'); - %% in case header has more than one row - header_str = string(join(append(join(header, ','), '\n'), '')); - + header = cellstr(header); % cellstr is for MATLAB - Octave compatibility + headerSize = size(header); + nHeaderLines = headerSize(1); % If there is a multi-line header + headerString = repmat({}, nHeaderLines, 1); + for k=1:nHeaderLines + % cellstr is for MATLAB - Octave compatibility + headerString(k) = cellstr(strjoin(header(k, :), ",")); + end + headerString = strcat(strjoin(headerString, "\n"), "\n"); + if not_header - header_str = [header_str]; + headerString = [headerString]; else % it is a header => each column must have one assert(length(header) == f_n_col, 'Less headers than lines `%s` or n_col is wrong', bin_path) end - fprintf(f_csv, header_str); + + fprintf(f_csv, headerString); fprintf(f_csv, [strjoin(units, ','), '\n']); - + f_bin = fopen(bin_path, 'r'); out = fread(f_bin, 'double'); -% fclose(f_bin); % + some useconds to execution - + out_2d = reshape(out, f_n_col, ns)'; -% dlmwrite(n_csv, out_2d, '-append', 'precision', '%d'); % SLOW! for k=1:ns fprintf(f_csv, '%d,', out_2d(k, 1:end-1)); fprintf(f_csv, '%d\n', out_2d(k, end)); % saves from extra comma end -% fclose(f_csv); end diff --git a/src/+io/create_output_files_binary.m b/src/+io/create_output_files_binary.m index a7936824..5b72afeb 100644 --- a/src/+io/create_output_files_binary.m +++ b/src/+io/create_output_files_binary.m @@ -1,114 +1,108 @@ -function [Output_dir, f, fnames] = create_output_files_binary(parameter_file, sitename, path_of_code, input_path, output_path, spectral, options) -%% Create Output dir -string = clock; -simulation_name = char(sitename); -outdir_name = sprintf('%s_%4.0f-%02.0f-%02.0f-%02.0f%02.0f', simulation_name, string(1:5)); -Output_dir = [fullfile(output_path, outdir_name) filesep]; +function [Output_dir, fnames] = create_output_files_binary(parameter_file, sitename, path_of_code, input_path, output_path, spectral, options) + %% Set Output dir + Output_dir = output_path; -warning('off','MATLAB:DELETE:FileNotFound') -if any(~exist(Output_dir,'dir')) - mkdir(Output_dir) + %% Log File mkdir([Output_dir,'Parameters' filesep]) -end + for i = 1:length(parameter_file) + copy_name = [strrep(parameter_file{i}, '.csv', '') '_' sitename '.csv']; + copyfile([input_path parameter_file{i}],[Output_dir,'Parameters/', copy_name],'f') + end + fidpath = fopen([Output_dir,'Parameters/SCOPEversion.txt'],'w'); % complete path of the SCOPE code + fprintf(fidpath,'%s', path_of_code); -%% Log File -for i = 1:length(parameter_file) - copy_name = [strrep(parameter_file{i}, '.csv', '') '_' outdir_name '.csv']; - copyfile([input_path parameter_file{i}],[Output_dir,'Parameters/', copy_name],'f') -end -fidpath = fopen([Output_dir,'Parameters/SCOPEversion.txt'],'w'); % complete path of the SCOPE code -fprintf(fidpath,'%s', path_of_code); + %% Filenames, will become .csv if options is on + fnames.flu_file = fullfile(Output_dir,'fluxes.bin'); + fnames.surftemp_file = fullfile(Output_dir,'surftemp.bin'); + fnames.aerodyn_file = fullfile(Output_dir,'aerodyn.bin'); + fnames.radiation_file = fullfile(Output_dir,'radiation.bin'); + fnames.fluorescence_file = fullfile(Output_dir,'fluorescence.bin'); + fnames.wl_file = fullfile(Output_dir,'wl.bin'); % wavelength + fnames.irradiance_spectra_file = fullfile(Output_dir,'irradiance_spectra.bin'); % Fluorescence + fnames.spectrum_hemis_optical_file = fullfile(Output_dir,'spectrum_hemis.bin'); + fnames.spectrum_obsdir_optical_file = fullfile(Output_dir,'spectrum_obsdir.bin'); + fnames.reflectance_file = fullfile(Output_dir,'reflectance.bin'); % reflectance spectrum + fnames.BOC_irradiance_file = fullfile(Output_dir,'BOC_irradiance.bin'); % reflectance spectrum + fnames.Sim_Theta_file = fullfile(Output_dir,'Sim_Theta.bin'); % soil moisture + fnames.Sim_Temp_file = fullfile(Output_dir,'Sim_Temp.bin'); % soil temperature + if options.calc_ebal + fnames.spectrum_obsdir_BlackBody_file = fullfile(Output_dir,'spectrum_obsdir_BlackBody.bin'); % spectrum observation direction + end -%% Filenames, will become .csv if options is on -fnames.flu_file = fullfile(Output_dir,'fluxes.bin'); -fnames.surftemp_file = fullfile(Output_dir,'surftemp.bin'); -fnames.aerodyn_file = fullfile(Output_dir,'aerodyn.bin'); -fnames.radiation_file = fullfile(Output_dir,'radiation.bin'); -fnames.fluorescence_file = fullfile(Output_dir,'fluorescence.bin'); -fnames.wl_file = fullfile(Output_dir,'wl.bin'); % wavelength -fnames.irradiance_spectra_file = fullfile(Output_dir,'irradiance_spectra.bin'); % Fluorescence -fnames.spectrum_hemis_optical_file = fullfile(Output_dir,'spectrum_hemis.bin'); -fnames.spectrum_obsdir_optical_file = fullfile(Output_dir,'spectrum_obsdir.bin'); -fnames.reflectance_file = fullfile(Output_dir,'reflectance.bin'); % reflectance spectrum -fnames.BOC_irradiance_file = fullfile(Output_dir,'BOC_irradiance.bin'); % reflectance spectrum -fnames.Sim_Theta_file = fullfile(Output_dir,'Sim_Theta.bin'); % soil moisture -fnames.Sim_Temp_file = fullfile(Output_dir,'Sim_Temp.bin'); % soil temperature -if options.calc_ebal - fnames.spectrum_obsdir_BlackBody_file = fullfile(Output_dir,'spectrum_obsdir_BlackBody.bin'); % spectrum observation direction -end + %if ~(options.simulation==1) + fnames.pars_and_input_file = fullfile(Output_dir,'pars_and_input.bin'); % wavelength + + %for j = 1:length(V) + % fprintf(fidv,'%s\t',V(j).Name); + %end + %fprintf(fidv,'\r'); + %end -%if ~(options.simulation==1) -fnames.pars_and_input_file = fullfile(Output_dir,'pars_and_input.bin'); % wavelength + %if ~(options.simulation==1) + fnames.pars_and_input_short_file = fullfile(Output_dir,'pars_and_input_short.bin'); % wavelength + %for j = find(vmax>1) + % fprintf(fidvs,'%s\t',V(vmax>1).Name); + %end + %fprintf(fidvs,' \r'); + % + %% Optional Output + if options.calc_vert_profiles + fnames.gap_file = fullfile(Output_dir,'gap.bin'); % gap + fnames.leaftemp_file = fullfile(Output_dir,'leaftemp.bin'); % leaftemp + fnames.layer_H_file = fullfile(Output_dir,'layer_H.bin'); % vertical profile + fnames.layer_LE_file = fullfile(Output_dir,'layer_lE.bin'); % latent heat + fnames.layer_A_file = fullfile(Output_dir,'layer_A.bin'); % + fnames.layer_aPAR_file = fullfile(Output_dir,'layer_aPAR.bin'); % + fnames.layer_aPAR_Cab_file = fullfile(Output_dir,'layer_aPAR_Cab.bin'); % + fnames.layer_Rn_file = fullfile(Output_dir,'layer_Rn.bin'); % + if options.calc_fluor + fnames.layer_fluorescence_file = fullfile(Output_dir,'layer_fluorescence.bin'); % + fnames.layer_fluorescenceEm_file = fullfile(Output_dir,'layer_fluorescenceEm.bin'); % + fnames.layer_NPQ_file = fullfile(Output_dir,'layer_NPQ.bin'); % + end -%for j = 1:length(V) - % fprintf(fidv,'%s\t',V(j).Name); -%end -%fprintf(fidv,'\r'); -%end + else + delete([Output_dir,'leaftemp.bin']) + delete([Output_dir,'layer_H.bin']) + delete([Output_dir,'layer_lE.bin']) + delete([Output_dir,'layer_A.bin']) + delete([Output_dir,'layer_aPAR.bin']) + delete([Output_dir,'layer_Rn.bin']) + end -%if ~(options.simulation==1) -fnames.pars_and_input_short_file = fullfile(Output_dir,'pars_and_input_short.bin'); % wavelength -%for j = find(vmax>1) - % fprintf(fidvs,'%s\t',V(vmax>1).Name); -%end -%fprintf(fidvs,' \r'); -% -%% Optional Output -if options.calc_vert_profiles - fnames.gap_file = fullfile(Output_dir,'gap.bin'); % gap - fnames.leaftemp_file = fullfile(Output_dir,'leaftemp.bin'); % leaftemp - fnames.layer_H_file = fullfile(Output_dir,'layer_H.bin'); % vertical profile - fnames.layer_LE_file = fullfile(Output_dir,'layer_lE.bin'); % latent heat - fnames.layer_A_file = fullfile(Output_dir,'layer_A.bin'); % - fnames.layer_aPAR_file = fullfile(Output_dir,'layer_aPAR.bin'); % - fnames.layer_aPAR_Cab_file = fullfile(Output_dir,'layer_aPAR_Cab.bin'); % - fnames.layer_Rn_file = fullfile(Output_dir,'layer_Rn.bin'); % if options.calc_fluor - fnames.layer_fluorescence_file = fullfile(Output_dir,'layer_fluorescence.bin'); % - fnames.layer_fluorescenceEm_file = fullfile(Output_dir,'layer_fluorescenceEm.bin'); % - fnames.layer_NPQ_file = fullfile(Output_dir,'layer_NPQ.bin'); % + fnames.fluorescence_file = fullfile(Output_dir,'fluorescence.bin'); % Fluorescence + if options.calc_PSI + fnames.fluorescencePSI_file = fullfile(Output_dir,'fluorescencePSI.bin'); % Fluorescence + fnames.fluorescencePSII_file = fullfile(Output_dir,'fluorescencePSII.bin'); % Fluorescence + end + fnames.fluorescence_hemis_file = fullfile(Output_dir,'fluorescence_hemis.bin'); % Fluorescence + fnames.fluorescence_emitted_by_all_leaves_file = fullfile(Output_dir,'fluorescence_emitted_by_all_leaves.bin'); + fnames.fluorescence_emitted_by_all_photosystems_file = fullfile(Output_dir,'fluorescence_emitted_by_all_photosystems.bin'); + fnames.fluorescence_sunlit_file = fullfile(Output_dir,'fluorescence_sunlit.bin'); % Fluorescence + fnames.fluorescence_shaded_file = fullfile(Output_dir,'fluorescence_shaded.bin'); % Fluorescence + fnames.fluorescence_scattered_file = fullfile(Output_dir,'fluorescence_scattered.bin'); % Fluorescence + else + delete([Output_dir,'fluorescence.bin']) end - -else - delete([Output_dir,'leaftemp.bin']) - delete([Output_dir,'layer_H.bin']) - delete([Output_dir,'layer_lE.bin']) - delete([Output_dir,'layer_A.bin']) - delete([Output_dir,'layer_aPAR.bin']) - delete([Output_dir,'layer_Rn.bin']) -end -if options.calc_fluor - fnames.fluorescence_file = fullfile(Output_dir,'fluorescence.bin'); % Fluorescence - if options.calc_PSI - fnames.fluorescencePSI_file = fullfile(Output_dir,'fluorescencePSI.bin'); % Fluorescence - fnames.fluorescencePSII_file = fullfile(Output_dir,'fluorescencePSII.bin'); % Fluorescence + if options.calc_directional + delete([Output_dir,'BRDF/*.bin']) end - fnames.fluorescence_hemis_file = fullfile(Output_dir,'fluorescence_hemis.bin'); % Fluorescence - fnames.fluorescence_emitted_by_all_leaves_file = fullfile(Output_dir,'fluorescence_emitted_by_all_leaves.bin'); - fnames.fluorescence_emitted_by_all_photosystems_file = fullfile(Output_dir,'fluorescence_emitted_by_all_photosystems.bin'); - fnames.fluorescence_sunlit_file = fullfile(Output_dir,'fluorescence_sunlit.bin'); % Fluorescence - fnames.fluorescence_shaded_file = fullfile(Output_dir,'fluorescence_shaded.bin'); % Fluorescence - fnames.fluorescence_scattered_file = fullfile(Output_dir,'fluorescence_scattered.bin'); % Fluorescence -else - delete([Output_dir,'fluorescence.bin']) -end -if options.calc_directional - delete([Output_dir,'BRDF/*.bin']) -end + if options.calc_planck && options.calc_ebal + fnames.spectrum_obsdir_thermal_file = fullfile(Output_dir,'spectrum_obsdir_thermal.bin'); % spectrum observation direction + fnames.spectrum_hemis_thermal_file = fullfile(Output_dir,'spectrum_hemis_thermal.bin'); % spectrum hemispherically integrated + end -if options.calc_planck && options.calc_ebal - fnames.spectrum_obsdir_thermal_file = fullfile(Output_dir,'spectrum_obsdir_thermal.bin'); % spectrum observation direction - fnames.spectrum_hemis_thermal_file = fullfile(Output_dir,'spectrum_hemis_thermal.bin'); % spectrum hemispherically integrated -end -%% Open files for writing -f = structfun(@(x) fopen(x, 'w'), fnames, 'UniformOutput',false); + % Create empty files for appending + structfun(@(x) fopen(x, 'w'), fnames, 'UniformOutput',false); + fclose("all"); -%% write wl -wlS = spectral.wlS; %#ok<*NASGU> -wlF = spectral.wlF; + %% write wl + wlS = spectral.wlS; %#ok<*NASGU> + wlF = spectral.wlF; -save([Output_dir 'wlS.txt'], 'wlS', '-ascii'); -save([Output_dir 'wlF.txt'], 'wlF', '-ascii'); + save([Output_dir 'wlS.txt'], 'wlS', '-ascii'); + save([Output_dir 'wlF.txt'], 'wlF', '-ascii'); end diff --git a/src/+io/output_data_binary.m b/src/+io/output_data_binary.m index 16dd97e0..0039b020 100644 --- a/src/+io/output_data_binary.m +++ b/src/+io/output_data_binary.m @@ -1,15 +1,7 @@ function n_col = output_data_binary(f, k, xyt, rad, canopy, V, vi, vmax, options, fluxes, meteo, iter, thermal,spectral, gap, profiles, Sim_Theta_U, Sim_Temp, Trap, Evap) %% OUTPUT DATA % author C. Van der Tol -% date: 30 Nov 2019 - -%% -if isdatetime(xyt.t) - get_doy = @(x) juliandate(x) - juliandate(datetime(year(x), 1, 0)); - V(46).Val = get_doy(timestamp2datetime(xyt.startDOY)); - V(47).Val = get_doy(timestamp2datetime(xyt.endDOY)); - xyt.t = get_doy(xyt.t); -end +% date: 30 Nov 2019 %% Fluxes product flu_out = [k iter.counter xyt.year(k) xyt.t(k) fluxes.Rntot fluxes.lEtot fluxes.Htot fluxes.Rnctot fluxes.lEctot, ... @@ -54,24 +46,24 @@ spectrum_obsdir_BlackBody_out = [rad.LotBB_']; n_col.spectrum_obsdir_BlackBody = length(spectrum_obsdir_BlackBody_out); fwrite(f.spectrum_obsdir_BlackBody_file,spectrum_obsdir_BlackBody_out,'double'); - + if options.calc_planck - - spectrum_hemis_thermal_out = [rad.Eoutte_']; + + spectrum_hemis_thermal_out = [rad.Eoutte_']; n_col.spectrum_hemis_thermal = length(spectrum_hemis_thermal_out); fwrite(f.spectrum_hemis_thermal_file,spectrum_hemis_thermal_out,'double'); - + spectrum_obsdir_thermal_out = [rad.Lot_']; n_col.spectrum_obsdir_thermal = length(spectrum_obsdir_thermal_out); fwrite(f.spectrum_obsdir_thermal_file,spectrum_obsdir_thermal_out,'double'); - + end end irradiance_spectra_out = [meteo.Rin*(rad.fEsuno+rad.fEskyo)']; n_col.irradiance_spectra = length(irradiance_spectra_out); fwrite(f.irradiance_spectra_file,irradiance_spectra_out,'double'); - + reflectance = pi*rad.Lo_./(rad.Esun_+rad.Esky_); reflectance(spectral.wlS>3000) = NaN; reflectance_out = [reflectance']; @@ -93,52 +85,52 @@ %% Optional Output if options.calc_vert_profiles - + % gap gap_out = [gap.Ps gap.Po gap.Pso]; n_col.gap = numel(gap_out(:)); fwrite(f.gap_file,gap_out,'double'); - + layer_aPAR_out = [1E6*profiles.Pn1d' 0]; n_col.layer_aPAR = length(layer_aPAR_out); fwrite(f.layer_aPAR_file,layer_aPAR_out,'double'); - + layer_aPAR_Cab_out = [1E6*profiles.Pn1d_Cab' 0]; n_col.layer_aPAR_Cab = length(layer_aPAR_Cab_out); - fwrite(f.layer_aPAR_Cab_file,layer_aPAR_Cab_out,'double'); - + fwrite(f.layer_aPAR_Cab_file,layer_aPAR_Cab_out,'double'); + if options.calc_ebal - - % leaftemp + + % leaftemp leaftemp_out = [profiles.Tcu1d' profiles.Tch' profiles.Tc1d']; n_col.leaftemp = length(leaftemp_out); - fwrite(f.leaftemp_file,leaftemp_out,'double'); - + fwrite(f.leaftemp_file,leaftemp_out,'double'); + layer_H_out = [profiles.Hc1d' fluxes.Hstot]; n_col.layer_H = length(layer_H_out); - fwrite(f.layer_H_file,layer_H_out,'double'); - + fwrite(f.layer_H_file,layer_H_out,'double'); + layer_LE_out = [profiles.lEc1d' fluxes.lEstot]; n_col.layer_LE = length(layer_LE_out); - fwrite(f.layer_LE_file,layer_LE_out,'double'); - + fwrite(f.layer_LE_file,layer_LE_out,'double'); + layer_A_out = [profiles.A1d' fluxes.Resp]; n_col.layer_A = length(layer_A_out); - fwrite(f.layer_A_file,layer_A_out,'double'); - + fwrite(f.layer_A_file,layer_A_out,'double'); + layer_NPQ_out = [profiles.qE' 0]; n_col.layer_NPQ = length(layer_NPQ_out); - fwrite(f.layer_NPQ_file,layer_NPQ_out,'double'); - + fwrite(f.layer_NPQ_file,layer_NPQ_out,'double'); + layer_Rn_out = [profiles.Rn1d' fluxes.Rnstot]; n_col.layer_Rn = length(layer_Rn_out); - fwrite(f.layer_Rn_file,layer_Rn_out,'double'); - + fwrite(f.layer_Rn_file,layer_Rn_out,'double'); + end if options.calc_fluor layer_fluorescence_out = [profiles.fluorescence']; n_col.layer_fluorescence = length(layer_fluorescence_out); - fwrite(f.layer_fluorescence_file,layer_fluorescence_out,'double'); + fwrite(f.layer_fluorescence_file,layer_fluorescence_out,'double'); end end @@ -146,47 +138,47 @@ for j=1:size(spectral.wlF,1) fluorescence_out = [rad.LoF_]; n_col.fluorescence = length(fluorescence_out); - fwrite(f.fluorescence_file,fluorescence_out,'double'); - + fwrite(f.fluorescence_file,fluorescence_out,'double'); + if options.calc_PSI fluorescencePSI_out = [rad.LoF1_]; n_col.fluorescencePSI = length(fluorescencePSI_out); - fwrite(f.fluorescencePSI_file,fluorescencePSI_out,'double'); - + fwrite(f.fluorescencePSI_file,fluorescencePSI_out,'double'); + fluorescencePSII_out = [rad.LoF2_]; n_col.fluorescencePSII = length(fluorescencePSII_out); - fwrite(f.fluorescencePSII_file,fluorescencePSII_out,'double'); + fwrite(f.fluorescencePSII_file,fluorescencePSII_out,'double'); end - + fluorescence_hemis_out = [rad.Fhem_]; n_col.fluorescence_hemis = length(fluorescence_hemis_out); - fwrite(f.fluorescence_hemis_file,fluorescence_hemis_out,'double'); - + fwrite(f.fluorescence_hemis_file,fluorescence_hemis_out,'double'); + fluorescence_emitted_by_all_leaves_out = [rad.Fem_]; n_col.fluorescence_emitted_by_all_leaves = length(fluorescence_emitted_by_all_leaves_out); - fwrite(f.fluorescence_emitted_by_all_leaves_file,fluorescence_emitted_by_all_leaves_out,'double'); - + fwrite(f.fluorescence_emitted_by_all_leaves_file,fluorescence_emitted_by_all_leaves_out,'double'); + fluorescence_emitted_by_all_photosystems_out = [rad.Femtot]; n_col.fluorescence_emitted_by_all_photosystems = length(fluorescence_emitted_by_all_photosystems_out); - fwrite(f.fluorescence_emitted_by_all_photosystems_file,fluorescence_emitted_by_all_photosystems_out,'double'); - + fwrite(f.fluorescence_emitted_by_all_photosystems_file,fluorescence_emitted_by_all_photosystems_out,'double'); + fluorescence_sunlit_out = [sum(rad.LoF_sunlit,2)]; n_col.fluorescence_sunlit = length(fluorescence_sunlit_out); - fwrite(f.fluorescence_sunlit_file,fluorescence_sunlit_out,'double'); - + fwrite(f.fluorescence_sunlit_file,fluorescence_sunlit_out,'double'); + fluorescence_shaded_out = [sum(rad.LoF_shaded,2)]; n_col.fluorescence_shaded = length(fluorescence_shaded_out); - fwrite(f.fluorescence_shaded_file,fluorescence_shaded_out,'double'); - + fwrite(f.fluorescence_shaded_file,fluorescence_shaded_out,'double'); + fluorescence_scattered_out = [sum(rad.LoF_scattered,2)+sum(rad.LoF_soil,2)]; n_col.fluorescence_scattered = length(fluorescence_scattered_out); - fwrite(f.fluorescence_scattered_file,fluorescence_scattered_out,'double'); - - end + fwrite(f.fluorescence_scattered_file,fluorescence_scattered_out,'double'); + + end end BOC_irradiance_out = [rad.Emin_(canopy.nlayers+1,:),rad.Emin_(canopy.nlayers+1,:)+(rad.Esun_*gap.Ps(canopy.nlayers+1)')']; n_col.BOC_irradiance = length(BOC_irradiance_out); - fwrite(f.BOC_irradiance_file,BOC_irradiance_out,'double'); + fwrite(f.BOC_irradiance_file,BOC_irradiance_out,'double'); %% if options.calc_directional && options.calc_ebal @@ -200,15 +192,15 @@ if options.calc_fluor Output_fluor = [spectral.wlF' directional.LoF_]; end - + save([Output_dir,'Directional/',sprintf('BRDF (SunAngle %2.2f degrees).dat',angles.tts)],'Output_brdf' ,'-ASCII','-TABS') save([Output_dir,'Directional/',sprintf('Angles (SunAngle %2.2f degrees).dat',angles.tts)],'Output_angle','-ASCII','-TABS') save([Output_dir,'Directional/',sprintf('Temperatures (SunAngle %2.2f degrees).dat',angles.tts)],'Output_temp','-ASCII','-TABS') - + if options.calc_fluor save([Output_dir,'Directional/',sprintf('Fluorescence (SunAngle %2.2f degrees).dat',angles.tts)],'Output_fluor','-ASCII','-TABS') end - + fiddirtir = fopen([Output_dir,'Directional/','read me.txt'],'w'); fprintf(fiddirtir,'The Directional data is written in three files: \r\n'); fprintf(fiddirtir,'\r\n- Angles: contains the directions. \r\n'); diff --git a/src/+io/prepareForcingData.m b/src/+io/prepareForcingData.m deleted file mode 100644 index e0f3c4a6..00000000 --- a/src/+io/prepareForcingData.m +++ /dev/null @@ -1,143 +0,0 @@ -function [SiteProperties, timeStep, forcingTimeLength] = prepareForcingData(DataPaths, forcingFileName) -%{ -This function is used to read forcing data and site properties. - -Input: - DataPaths: A construct contains paths of data. - forcingFileName: A string of the name of forcing NetCDF data. - -Output: - SiteProperties： A structure containing site properties variables. - timeStep: Time interval in seconds, normal is 1800 s in STEMMUS-SCOPE. - timeLength: The total number of time steps in forcing data. -%} -%% -% Add the forcing name to forcing path -ForcingFilePath=fullfile(DataPaths.forcingPath, forcingFileName); -% Prepare input files -sitefullname=dir(ForcingFilePath).name; %read sitename -SiteProperties.siteName=sitefullname(1:6); -startyear=sitefullname(8:11); -endyear=sitefullname(13:16); -startyear=str2double(startyear); -endyear=str2double(endyear); - -% Read time values from forcing file -time1=ncread(ForcingFilePath,'time'); -t1=datenum(startyear,1,1,0,0,0); -timeStep=time1(2); - -%get time length of forcing file -forcingTimeLength=length(time1); - -dt=time1(2)/3600/24; -t2=datenum(endyear,12,31,23,30,0); -T=t1:dt:t2; -TL=length(T); -T=T'; -T=datestr(T,'yyyy-mm-dd HH:MM:SS'); -T0=T(:,1:4); -T1=T(:,5:19); -T3=T1(1,:); -T4=T3(ones(TL,1),:); -T5=[T0,T4]; -T6=datenum(T); -T7=datenum(T5); -T8=T6-T7; -time=T8; -T10=year(T); - -RH=ncread(ForcingFilePath,'RH'); % Unit: % -RHL=length(RH); -RHa=reshape(RH,RHL,1); - -Tair=ncread(ForcingFilePath,'Tair'); % Unit: K -TairL=length(Tair); -Taira=reshape(Tair,TairL,1)-273.15; % Unit: degree C - -es= 6.107*10.^(Taira.*7.5./(237.3+Taira)); % Unit: hPa -ea=es.*RHa./100; - -SWdown=ncread(ForcingFilePath,'SWdown'); % Unit: W/m2 -SWdownL=length(SWdown); -SWdowna=reshape(SWdown,SWdownL,1); - - -LWdown=ncread(ForcingFilePath,'LWdown'); % Unit: W/m2 -LWdownL=length(LWdown); -LWdowna=reshape(LWdown,LWdownL,1); - - -VPD=ncread(ForcingFilePath,'VPD'); % Unit: hPa -VPDL=length(VPD); -VPDa=reshape(VPD,VPDL,1); - - -% Qair=ncread(ForcingFilePath,'Qair'); -% QairL=length(Qair); -% Qaira=reshape(Qair,QairL,1); - - -Psurf=ncread(ForcingFilePath,'Psurf'); % Unit: Pa -PsurfL=length(Psurf); -Psurfa=reshape(Psurf,PsurfL,1); -Psurfa=Psurfa./100; % Unit: hPa - - -Precip=ncread(ForcingFilePath,'Precip'); % Unit: mm/s -PrecipL=length(Precip); -Precipa=reshape(Precip,PrecipL,1); -Precipa=Precipa./10; % Unit: cm/s - - -Wind=ncread(ForcingFilePath,'Wind'); % Unit: m/s -WindL=length(Wind); -Winda=reshape(Wind,WindL,1); - - -CO2air=ncread(ForcingFilePath,'CO2air'); % Unit: ppm -CO2airL=length(CO2air); -CO2aira=reshape(CO2air,CO2airL,1); -CO2aira=CO2aira.*44./22.4; % Unit: mg/m3 - - -SiteProperties.latitude=ncread(ForcingFilePath,'latitude'); -SiteProperties.longitude=ncread(ForcingFilePath,'longitude'); -SiteProperties.elevation=ncread(ForcingFilePath,'elevation'); - -LAI=ncread(ForcingFilePath,'LAI'); -LAIL=length(LAI); -LAIa=reshape(LAI,LAIL,1); -LAIa(LAIa<0.01)=0.01; - -% LAI_alternative=ncread(ForcingFilePath,'LAI_alternative'); -% LAI_alternativeL=length(LAI_alternative); -% LAI_alternativea=reshape(LAI_alternative,LAI_alternativeL,1); - -SiteProperties.igbpVegLong=ncread(ForcingFilePath,'IGBP_veg_long'); -SiteProperties.referenceHeight=ncread(ForcingFilePath,'reference_height'); -SiteProperties.canopyHeight=ncread(ForcingFilePath,'canopy_height'); - -save([DataPaths.input, 't_.dat'], '-ascii', 'time') -save([DataPaths.input, 'Ta_.dat'], '-ascii', 'Taira') -save([DataPaths.input, 'Rin_.dat'], '-ascii', 'SWdowna') -save([DataPaths.input, 'Rli_.dat'], '-ascii', 'LWdowna') -%save([DataPaths.input, 'VPDa.dat'], '-ascii', 'VPDa') -%save([DataPaths.input, 'Qaira.dat'], '-ascii', 'Qaira') -save([DataPaths.input, 'p_.dat'], '-ascii', 'Psurfa') -%save([DataPaths.input, 'Precipa.dat'], '-ascii', 'Precipa') -save([DataPaths.input, 'u_.dat'], '-ascii', 'Winda') -%save([DataPaths.input, 'RHa.dat'], '-ascii', 'RHa') -save([DataPaths.input, 'CO2_.dat'], '-ascii', 'CO2aira') -%save([DataPaths.input, 'latitude.dat'], '-ascii', 'latitude') -%save([DataPaths.input, 'longitude.dat'], '-ascii', 'longitude') -%save([DataPaths.input, 'reference_height.dat'], '-ascii', 'reference_height') -%save([DataPaths.input, 'canopy_height.dat'], '-ascii', 'canopy_height') -%save([DataPaths.input, 'elevation.dat'], '-ascii', 'elevation') -save([DataPaths.input, 'ea_.dat'], '-ascii', 'ea') -save([DataPaths.input, 'year_.dat'], '-ascii', 'T10') -LAI=[time'; LAIa']'; -save([DataPaths.input, 'LAI_.dat'], '-ascii', 'LAI') %save meteorological data for STEMMUS -Meteodata=[time';Taira';RHa';Winda';Psurfa';Precipa';SWdowna';LWdowna';VPDa';LAIa']'; -save([DataPaths.input, 'Mdata.txt'], '-ascii', 'Meteodata') %save meteorological data for STEMMUS -end \ No newline at end of file diff --git a/src/+io/prepareInputData.m b/src/+io/prepareInputData.m new file mode 100644 index 00000000..11b2bec4 --- /dev/null +++ b/src/+io/prepareInputData.m @@ -0,0 +1,34 @@ +function [SiteProperties, SoilProperties, TimeProperties] = prepareInputData(InputPath) + %{ + This function is used to read forcing data and site properties. + + Input: + InputPath: path of input data. + + Output: + SiteProperties： A structure containing site properties variables. + SoilProperties: A structure containing soil variables. + TimeProperties: A structure containing time variables like time interval in seconds, normal is 1800 s in STEMMUS-SCOPE + and the total number of time steps. + %} + + % The "forcing_globals.mat" and "soil_parameters.mat" files are generated using "PyStemmusScope" + % python package, see https://github.com/EcoExtreML/STEMMUS_SCOPE_Processing + forcing_global_path = fullfile(InputPath, 'forcing_globals.mat'); + soil_global_path = fullfile(InputPath, 'soil_parameters.mat'); + + %% Explicitly load all variables into the workspace + TimeProperties = load(forcing_global_path, 'DELT', 'Dur_tot'); + SiteProperties = load(forcing_global_path, 'IGBP_veg_long', 'latitude', 'longitude', 'reference_height', 'canopy_height', 'sitename'); + + SoilProperties = load(soil_global_path, 'SaturatedK', 'SaturatedMC', 'ResidualMC', 'Coefficient_n', 'Coefficient_Alpha', ... + 'porosity', 'FOC', 'FOS', 'MSOC', 'Coef_Lamda', 'fieldMC', 'fmax', 'theta_s0', 'Ks0'); + + % Convert the int vectors back to strings + SiteProperties.sitename = char(SiteProperties.sitename); + SiteProperties.IGBP_veg_long = char(SiteProperties.IGBP_veg_long); + + % The model expects the char vector to be transposed: + SiteProperties.IGBP_veg_long = transpose(SiteProperties.IGBP_veg_long); + +end \ No newline at end of file diff --git a/src/+io/readStructFromExcel.m b/src/+io/readStructFromExcel.m index 30e0ff2a..210ac7a8 100644 --- a/src/+io/readStructFromExcel.m +++ b/src/+io/readStructFromExcel.m @@ -1,52 +1,53 @@ function data = readStructFromExcel(filename, sheetName, headerIdx, dataIdx, data_is_char, data_in_rows) -% Read data into a struct with names matching those found in the first column/row -% default is for data to be in columns (A and B); if data_in_rows = true, data are in rows 1 & 2 -% example: -% readStructFromExcel('../input_data.xlsx', 'options', 3, 1) -% readStructFromExcel('../input_data.xlsx', 'filenames', 1, 2, true) + % Read data into a struct with names matching those found in the first column/row + % default is for data to be in columns (A and B); if data_in_rows = true, data are in rows 1 & 2 + % example: + % readStructFromExcel('../input_data.xlsx', 'options', 3, 1) + % readStructFromExcel('../input_data.xlsx', 'filenames', 1, 2, true) -% default values: -if nargin < 3 - headerIdx = 1; -end -if nargin < 4 - dataIdx = 2; -end -if nargin < 5 - data_is_char = false; -end -if nargin < 6 - data_in_rows = false; -end + % default values: + if nargin < 3 + headerIdx = 1; + end + if nargin < 4 + dataIdx = 2; + end + if nargin < 5 + data_is_char = false; + end + if nargin < 6 + data_in_rows = false; + end -% read in the spreadsheet -% general note: work with all_as_cell to keep rows & columns in sync; MATLAB does NOT keep data and texts aligned -% (readtable only slightly better - it works but will convert mixed columns to string) -if data_in_rows - %NOTE: 'basic' is compatible with Mac but MATLAB (2013b) complains it can't read Unicode '.xls' in basic mode - % Solution: save as xlsx ?! - [~, ~, allCells] = xlsread(filename, sheetName, ''); % = [data, texts, allCells] -else - [~, ~, allCells] = xlsread(filename, sheetName, ''); % = [data, texts, allCells] - allCells = allCells'; % transpose so data are in the same column as headers -end -% data are now in columns + % read in the spreadsheet + % general note: work with all_as_cell to keep rows & columns in sync; MATLAB does NOT keep data and texts aligned + % (readtable only slightly better - it works but will convert mixed columns to string) + if data_in_rows + %NOTE: 'basic' is compatible with Mac but MATLAB (2013b) complains it can't read Unicode '.xls' in basic mode + % Solution: save as xlsx ?! + [~, ~, allCells] = xlsread(filename, sheetName); % = [data, texts, allCells] + else + [~, ~, allCells] = xlsread(filename, sheetName); % = [data, texts, allCells] + allCells = allCells'; % transpose so data are in the same column as headers + end + % data are now in columns -% delete empty columns - % define two "helper functions" for eliminating null entries -isNotNumeric = @(x) cellfun(@(y) ischar(y) | any(isnan(y)), x); % any is needed because matlab treats string as char array -isCharCell = @(x) cellfun(@(y) ischar(y), x); + % delete empty columns + % define two "helper functions" for eliminating null entries + isNotNumeric = @(x) cellfun(@(y) ischar(y) | any(isnan(y)), x); % any is needed because matlab treats string as char array + isCharCell = @(x) cellfun(@(y) ischar(y), x); -validHeaders = arrayfun(isCharCell, allCells(headerIdx, :)); -if data_is_char - validData = arrayfun(isCharCell, allCells(dataIdx, :)); % , 'UniformOutput', true -else % numeric data: - validData = ~arrayfun(isNotNumeric, allCells(dataIdx, :)); % , 'UniformOutput', true -end + validHeaders = arrayfun(isCharCell, allCells(headerIdx, :)); + if data_is_char + validData = arrayfun(isCharCell, allCells(dataIdx, :)); % , 'UniformOutput', true + else % numeric data: + validData = ~arrayfun(isNotNumeric, allCells(dataIdx, :)); % , 'UniformOutput', true + end -dataCells = allCells([headerIdx, dataIdx], validHeaders & validData); + dataCells = allCells([headerIdx, dataIdx], validHeaders & validData); -for idx = 1:size(dataCells,2) - varName = strrep(dataCells{1, idx}, ' ', ''); - data.(varName) = dataCells{2, idx}; + for idx = 1:size(dataCells,2) + varName = strrep(dataCells{1, idx}, ' ', ''); + data.(varName) = dataCells{2, idx}; + end end \ No newline at end of file diff --git a/src/+io/read_config.m b/src/+io/read_config.m index a8f3840c..acf1eef5 100644 --- a/src/+io/read_config.m +++ b/src/+io/read_config.m @@ -1,4 +1,4 @@ -function [SoilPropertyPath, InputPath, OutputPath, ForcingPath, ForcingFileName, DurationSize, InitialConditionPath] = read_config(config_file) +function [InputPath, OutputPath, InitialConditionPath] = read_config(config_file) file_id = fopen(config_file); config = textscan(file_id,'%s %s', 'HeaderLines',0, 'Delimiter', '='); @@ -9,24 +9,11 @@ config_paths = config{2}; %% find the required path by model -indx = find(strcmp(config_vars, 'SoilPropertyPath')); -SoilPropertyPath = config_paths{indx}; - indx = find(strcmp(config_vars, 'InputPath')); InputPath = config_paths{indx}; indx = find(strcmp(config_vars, 'OutputPath')); OutputPath = config_paths{indx}; -indx = find(strcmp(config_vars, 'ForcingPath')); -ForcingPath = config_paths{indx}; - -indx = find(strcmp(config_vars, 'ForcingFileName')); -ForcingFileName = config_paths{indx}; - indx = find(strcmp(config_vars, 'InitialConditionPath')); InitialConditionPath = config_paths{indx}; - -% value of DurationSize is optional and can be NA -indx = find(strcmp(config_vars, 'DurationSize')); -DurationSize = str2double(config_paths{indx}); diff --git a/src/Constants.m b/src/Constants.m index d14fd0f2..194d1503 100644 --- a/src/Constants.m +++ b/src/Constants.m @@ -1,73 +1,73 @@ % function Constants global DeltZ Delt_t ML NS mL mN nD NL NN SAVE Tot_Depth Tmin global xERR hERR TERR PERR tS uERR -global KT TIME Delt_t0 DURTN TEND NIT KIT Nmsrmn Eqlspace h_SUR Msrmn_Fitting +global KT TIME Delt_t0 DURTN TEND NIT KIT Nmsrmn Eqlspace h_SUR Msrmn_Fitting global Msr_Mois Msr_Temp Msr_Time Tp_t Evap EPCT SAVEDTheta_LLT SAVEDTheta_LLh SAVEDTheta_UUh -global Ksoil Rl SMC Ztot DeltZ_R Theta_o rroot frac bbx wfrac RWUtot RWUtott RWUtottt Rls Tatot LR PSItot sfactortot LAI_msr R_depth RTB VPD_msr Tss Tsss SUMTIME Tsur FOC FOS FOSL MSOC Coef_Lamda fieldMC DELT Dur_tot Precipp fmax sitename +global Ksoil Rl SMC Ztot DeltZ_R Theta_o rroot frac bbx wfrac RWUtot RWUtott RWUtottt Rls Tatot LR PSItot sfactortot LAI_msr R_depth RTB VPD_msr Tss Tsss SUMTIME Tsur Precipp %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%% The time and domain information setting. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -KIT=0; % KIT is used to count the number of iteration in a time step; -NIT=30; % Desirable number of iterations in a time step; -Nmsrmn=Dur_tot*10; %Nmsrmn=140256*100; Here, it is made as big as possible, in case a long simulation period containing many time step is defined. -SUMTIME=0; -DURTN=DELT*Dur_tot; % Duration of simulation period; -KT=0; % Number of time steps; -TIME=0*DELT; % Time of simulation released; +KIT=0; % KIT is used to count the number of iteration in a time step; +NIT=30; % Desirable number of iterations in a time step; +Nmsrmn=Dur_tot*10; %Nmsrmn=140256*100; Here, it is made as big as possible, in case a long simulation period containing many time step is defined. +SUMTIME=0; +DURTN=DELT*Dur_tot; % Duration of simulation period; +KT=0; % Number of time steps; +TIME=0*DELT; % Time of simulation released; - -TEND=TIME+DURTN; % Time to be reached at the end of simulation period; -Delt_t=DELT; % Duration of time step [Unit of second] -Delt_t0=Delt_t; % Duration of last time step; -tS=DURTN/Delt_t; % Is the tS(time step) needed to be added with 1? + +TEND=TIME+DURTN; % Time to be reached at the end of simulation period; +Delt_t=DELT; % Duration of time step [Unit of second] +Delt_t0=Delt_t; % Duration of last time step; +tS=DURTN/Delt_t; % Is the tS(time step) needed to be added with 1? % Cause the start of simulation period is from 0mins, while the input data start from 30mins. - -xERR=0.02; % Maximum desirable change of moisture content; -hERR=0.1e08; % Maximum desirable change of matric potential; -TERR=2; % Maximum desirable change of temperature; + +xERR=0.02; % Maximum desirable change of moisture content; +hERR=0.1e08; % Maximum desirable change of matric potential; +TERR=2; % Maximum desirable change of temperature; PERR=5000; % Maximum desirable change of soil air pressure (Pa,kg.m^-1.s^-1); -uERR=0.02; % Maximum desirable change of total water content; -Tot_Depth=500; % Unit is cm. it should be usually bigger than 0.5m. Otherwise, - % the DeltZ would be reset in 50cm by hand; -R_depth=300; % -Eqlspace=0; % Indicator for deciding is the space step equal or not; +uERR=0.02; % Maximum desirable change of total water content; +Tot_Depth=500; % Unit is cm. it should be usually bigger than 0.5m. Otherwise, + % the DeltZ would be reset in 50cm by hand; +R_depth=300; % +Eqlspace=0; % Indicator for deciding is the space step equal or not; NL=100; -if ~Eqlspace - run Dtrmn_Z % Determination of NL, the number of elments; +if ~Eqlspace + run Dtrmn_Z % Determination of NL, the number of elments; else - for ML=1:NL - DeltZ(ML)=Tot_Depth/NL; - end -end - -NN=NL+1; % Number of nodes; -mN=NN+1; + for ML=1:NL + DeltZ(ML)=Tot_Depth/NL; + end +end + +NN=NL+1; % Number of nodes; +mN=NN+1; mL=NL+1; % Number of elements. Prevending the exceeds of size of arraies; -nD=2; -NS=1; % Number of soil types; -SAVE=zeros(3,3,3); % Arraies for calculating boundary flux; +nD=2; +NS=1; % Number of soil types; +SAVE=zeros(3,3,3); % Arraies for calculating boundary flux; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -global h_TE -global W_Chg +global h_TE +global W_Chg global l CKTN -global RHOL SSUR RHO_bulk Rv g RDA -global KL_h KL_T D_Ta +global RHOL SSUR RHO_bulk Rv g RDA +global KL_h KL_T D_Ta global Theta_L Theta_LL Se h hh T TT Theta_V h_frez hh_frez T_CRIT TT_CRIT global W WW MU_W f0 L_WT -global DhT +global DhT global GWT Gamma0 MU_W0 MU1 b W0 Gamma_w global Chh ChT Khh KhT Kha Vvh VvT Chg global C1 C2 C3 C4 C5 C6 C7 C9 global QL QL_D QL_disp RHS global HR RHOV_s RHOV DRHOV_sT DRHOVh DRHOVT global RHODA DRHODAt DRHODAz Xaa XaT Xah -global D_Vg D_V D_A +global D_Vg D_V D_A global k_g Sa V_A Alpha_Lg POR_C Eta global P_g P_gg Theta_g P_g0 Beta_g global MU_a fc Unit_C Hc UnitC global Cah CaT Caa Kah KaT Kaa Vah VaT Vaa Cag global Lambda1 Lambda2 Lambda3 c_L c_a c_V c_i L0 -global Lambda_eff c_unsat L LL Tr +global Lambda_eff c_unsat L LL Tr global CTh CTT CTa KTh KTT KTa VTT VTh VTa CTg CTT_PH CTT_LT CTT_g CTT_Lg EfTCON TETCON global Kcvh KcvT Kcva Ccvh CcvT Kcah KcaT Kcaa Ccah CcaT Ccaa @@ -76,14 +76,14 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% global MO Ta U Ts Zeta_MO % U_wind is the mean wind speed at height z_ref (m��s^-1), U is the wind speed at each time step. global Precip SH HR_a UseTs_msrmn % Notice that Evap and Precip have only one value for each time step. Precip needs to be measured in advance as the input. -global Gsc Sigma_E -global Rns Rnl -global Tmax Tmin Hrmax Hrmin Um +global Gsc Sigma_E +global Rns Rnl +global Tmax Tmin Hrmax Hrmin Um %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Variables information for soil air pressure propagation +% Variables information for soil air pressure propagation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -global TopPg +global TopPg %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Fluxes information with different mechanisms @@ -94,8 +94,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Variable information for updating the state variables %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -global hOLD TOLD P_gOLD J -global porosity SaturatedMC ResidualMC Coefficient_n Coefficient_Alpha +global hOLD TOLD P_gOLD J %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Variables information for initialization subroutine %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -108,22 +107,22 @@ % Effective Thermal conductivity and capacity variables %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% global EHCAP ThmrlCondCap Evaptranp_Cal SWCC ThermCond -global Ts_a0 Ts_a Ts_b Ts_w -global Pg_w Pg_a0 Pg_a Pg_b +global Ts_a0 Ts_a Ts_b Ts_w +global Pg_w Pg_a0 Pg_a Pg_b global Hsur_w Hsur_a0 Hsur_a Hsur_b Rn VPD_msr LAI_msr %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Variables information for boundary condition settings %> NBCh Indicator for type of surface boundary condition on mass euqation to be applied; -%> "1"--Specified matric head; +%> "1"--Specified matric head; %> "2"--Specified potential moisture flux (potential evaporation rate of precipitation rate); %> "3"--Specified atmospheric forcing; %> NBCT Indicator for type of surface boundary condtion on energy equation to be applied; -%> "1"--Specified temperature; +%> "1"--Specified temperature; %> "2"--Specified heat flux; % "3"--Specified atmospheric forcing; %> NBCP Indicator for type of surface boundary condition on dry air equation to be applied; -%> "1"--Specified atmospheric pressure; -%> "2"--Specified atmospheric forcing (Measured atmospheric pressure data); +%> "1"--Specified atmospheric pressure; +%> "2"--Specified atmospheric forcing (Measured atmospheric pressure data); %> BCh Value of specified matric head (if NBCh=1) or specified potential moisture flux (if NBCh=2) at surface; %> BChB The same as BCh but at the bottom of column; %> BCT Value of specified temperature (if NBCT=1) or specified head flux (if NBCT=2) at surface; @@ -136,7 +135,7 @@ %> "2"--Specified moisture flux of BChB; %> "3"--Zero matric head gradient at bottom(Gravity drainage); %> NBCTB Type of boundary condition on temperature at bottom of column; -%> "1"--Specified temperature BCTB at bottom; +%> "1"--Specified temperature BCTB at bottom; %> "2"--Specified heat flux BCTB at bottom; %> "3"--Zero temperature gradient at bottom (advection only); %> NBCPB Type of boundary condition on soil air pressure; @@ -150,7 +149,7 @@ %> DSTMAX Depression storage capacity; %> EXCESS Excess of available water over infiltration rate for current time step,expressed as a rate; %> AVAIL Maximum rate at which water can be supplied to the soil from above during the current time step; -%> +%> %> global alpha_h bx Srt rwuef Ts_msr Tb_msr Theta_U Theta_UU Theta_I Theta_II KfL_h RHOI KfL_T Tbtm r_a_SOIL Rn_SOIL EVAP XCAP SAVEKfL_h global Ta_msr RH_msr Rn_msr WS_msr G_msr Pg_msr HourInput Rns_msr SH_msr LE_msr DRHOVZ DTDZ DhDZ Precip_msr LAI_msr DTheta_LLh SAVETheta_UU DTheta_UUh EVAPO %SAVEDTheta_UUh @@ -158,34 +157,34 @@ bx=zeros(mL,nD);SAVEKfL_h=zeros(mL,nD); Srt=zeros(mL,nD);DhDZ=zeros(mL,1);DTDZ=zeros(mL,1);DRHOVZ=zeros(mL,1); DTheta_LLh=zeros(mL,nD);DTheta_UUh=zeros(mL,nD);SAVETheta_UU=zeros(mL,nD);SAVEDTheta_UUh=zeros(mL,nD);SAVEDTheta_LLh=zeros(mL,nD);SAVEDTheta_LLT=zeros(mL,nD); Ratio_ice=zeros(mL,nD); -KL_h=zeros(mL,nD); % The hydraulic conductivity(m��s^-1); -KfL_h=zeros(mL,nD); % The hydraulic conductivity considering ice blockking effect(m��s^-1); -KfL_T=zeros(mL,nD); % The depression temperature controlled by ice(m^2��Cels^-1��s^-1); -KL_T=zeros(mL,nD); % The conductivity controlled by thermal gradient(m^2��Cels^-1��s^-1); -D_Ta=zeros(mL,nD); % The thermal dispersivity for soil water (m^2��Cels^-1��s^-1); -Theta_L=zeros(mL,nD); % The soil moisture at the start of current time step; -Theta_LL=zeros(mL,nD); % The soil moisture at the end of current time step; -Theta_U=zeros(mL,nD); % The total soil moisture(water+ice) at the start of current time step; -Theta_UU=zeros(mL,nD); % The total soil moisture at the end of current time step; -Theta_I=zeros(mL,nD); % The soil ice water content at the start of current time step; -Theta_II=zeros(mL,nD); % The soil ice water content at the end of current time step; -Se=zeros(mL,nD); % The saturation degree of soil moisture; +KL_h=zeros(mL,nD); % The hydraulic conductivity(m��s^-1); +KfL_h=zeros(mL,nD); % The hydraulic conductivity considering ice blockking effect(m��s^-1); +KfL_T=zeros(mL,nD); % The depression temperature controlled by ice(m^2��Cels^-1��s^-1); +KL_T=zeros(mL,nD); % The conductivity controlled by thermal gradient(m^2��Cels^-1��s^-1); +D_Ta=zeros(mL,nD); % The thermal dispersivity for soil water (m^2��Cels^-1��s^-1); +Theta_L=zeros(mL,nD); % The soil moisture at the start of current time step; +Theta_LL=zeros(mL,nD); % The soil moisture at the end of current time step; +Theta_U=zeros(mL,nD); % The total soil moisture(water+ice) at the start of current time step; +Theta_UU=zeros(mL,nD); % The total soil moisture at the end of current time step; +Theta_I=zeros(mL,nD); % The soil ice water content at the start of current time step; +Theta_II=zeros(mL,nD); % The soil ice water content at the end of current time step; +Se=zeros(mL,nD); % The saturation degree of soil moisture; h=zeros(mN,1); % The matric head at the start of current time step; hh=zeros(mN,1); % The matric head at the end of current time step; h_frez=zeros(mN,1); % The freeze depression head at the start of current time step; hh_frez=zeros(mN,1); % The freeze depression head at the end of current time step; XCAP=zeros(mN,1); -T=zeros(mN,1); % The soil temperature at the start of current time step; -TT=zeros(mN,1); % The soil temperature at the end of current time step; -T_CRIT=zeros(mN,1); % The soil ice critical temperature at the start of current time step; -TT_CRIT=zeros(mN,1); % The soil ice critical temperature at the start of current time step; +T=zeros(mN,1); % The soil temperature at the start of current time step; +TT=zeros(mN,1); % The soil temperature at the end of current time step; +T_CRIT=zeros(mN,1); % The soil ice critical temperature at the start of current time step; +TT_CRIT=zeros(mN,1); % The soil ice critical temperature at the start of current time step; EPCT=zeros(mN,1); -Theta_V=zeros(mL,nD); % Volumetric gas content; -W=zeros(mL,nD); % Differential heat of wetting at the start of current time step(J��kg^-1); +Theta_V=zeros(mL,nD); % Volumetric gas content; +W=zeros(mL,nD); % Differential heat of wetting at the start of current time step(J��kg^-1); WW=zeros(mL,nD); % Differential heat of wetting at the end of current time step(J��kg^-1); % Integral heat of wetting in individual time step(J��m^-2); %%%%%%%%%%%%%%% Notice: the formulation of this in 'CondL_Tdisp' is not a sure. %%%%%%%%%%%%%% MU_W=zeros(mL,nD); % Visocity of water(kg��m^?6?1��s^?6?1); -f0=zeros(mL,nD); % Tortusity factor [Millington and Quirk (1961)]; kg.m^2.s^-2.m^-2.kg.m^-3 +f0=zeros(mL,nD); % Tortusity factor [Millington and Quirk (1961)]; kg.m^2.s^-2.m^-2.kg.m^-3 L_WT=zeros(mL,nD); % Liquid dispersion factor in Thermal dispersivity(kg��m^-1��s^-1)=-------------------------- m^2 (1.5548e-013 m^2); DhT=zeros(mN,1); % Difference of matric head with respect to temperature; m. kg.m^-1.s^-1 RHS=zeros(mN,1); % The right hand side part of equations in '*_EQ' subroutine; @@ -252,8 +251,8 @@ CTT_Lg=zeros(mL,nD); % Storage coefficient in energy conservation equation related to liquid and gas; CTT_g=zeros(mL,nD); % Storage coefficient in energy conservation equation related to air; CTT_LT=zeros(mL,nD); % Storage coefficient in energy conservation equation related to liquid and temperature; -EfTCON=zeros(mL,nD); % Effective heat conductivity Johansen method; -TETCON=zeros(mL,nD); % Effective heat conductivity Johansen method; +EfTCON=zeros(mL,nD); % Effective heat conductivity Johansen method; +TETCON=zeros(mL,nD); % Effective heat conductivity Johansen method; L=zeros(mN,1); % The latent heat of vaporization at the beginning of the time step; LL=zeros(mN,1); % The latent heat of vaporization at the end of the time step; CTh=zeros(mL,nD); % Storage coefficient in energy conservation equation related to matric head; @@ -325,7 +324,7 @@ % DTDZ=zeros(mL,1); KLhBAR=zeros(mL,1); DEhBAR=zeros(mL,1); -KLTBAR=zeros(mL,1); +KLTBAR=zeros(mL,1); DTDBAR=zeros(mL,1); QLH=zeros(mL,1); QLT=zeros(mL,1); @@ -357,34 +356,34 @@ Tsur=zeros(Nmsrmn/10,1); Ztot=zeros(ML,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%%%%%%% The indicators needs to be set before the running of this program %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%% The indicators needs to be set before the running of this program %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% J=1; % Indicator denotes the index of soil type for choosing soil physical parameters; % rwuef=1; HourInput=1; SWCC=1; % indicator for choose the soil water characteristic curve, =0, Clapp and Hornberger; =1, Van Gen; -Evaptranp_Cal=3; % Indicator denotes the method of estimating evapotranspiration; - % Value of 1 means the ETind method, otherwise, ETdir method; +Evaptranp_Cal=3; % Indicator denotes the method of estimating evapotranspiration; + % Value of 1 means the ETind method, otherwise, ETdir method; UseTs_msrmn=1; % Value of 1 means the measurement Ts would be used; Otherwise, 0; % -Msrmn_Fitting=1; % Value of 1 means the measurement data is used to fit the simulations; +Msrmn_Fitting=1; % Value of 1 means the measurement data is used to fit the simulations; Hystrs=0; % If the value of Hystrs is 1, then the hysteresis is considered, otherwise 0; % Thmrlefc=1; % Consider the isothermal water flow if the value is 0, otherwise 1; % Soilairefc=0; % The dry air transport is considered with the value of 1,otherwise 0; % hThmrl=1; % Value of 1, the special calculation of water capacity is used, otherwise 0; % h_TE=0; % Value of 1 means that the temperature dependence % - % of matric head would be considered.Otherwise,0; % + % of matric head would be considered.Otherwise,0; % W_Chg=1; % Value of 0 means that the heat of wetting would % - % be calculated by Milly's method��Otherwise,1. The % + % be calculated by Milly's method��Otherwise,1. The % % method of Lyle Prunty would be used; % ThmrlCondCap=1; %1; % The indicator for choosing Milly's effective thermal capacity and conductivity % - % formulation to verify the vapor and heat transport in extremly dry soil. % + % formulation to verify the vapor and heat transport in extremly dry soil. % % ThmrlCond=1; % The indicator for choosing effective thermal conductivity methods, 1= de vries method;2= Jonhansen methods ThermCond=1; % The indicator for choosing effective thermal conductivity methods, 1= de vries method;2= Jonhansen methods;3= Simplified de vries method(Tian 2016);4= Farouki methods % SWCC=1; % indicator for choose the soil water characteristic curve, =0, Clapp and Hornberger; =1, Van Gen; CHST=0; % Indicator of parameters derivation using soil texture or not. CHST=1, use; CHST=0 not use ISOC=1; % Indicator of considering soil organic matter effect or not. ISOC=1, yes; ISOC=0 no %%%%% 172 27%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -run soilpropertyread %load soil property + CKTN=(50+2.575*20); % Constant used in calculating viscosity factor for hydraulic conductivity l=0.5; % Coefficient in VG model; g=981; % Gravity acceleration (cm.s^-2); @@ -399,9 +398,9 @@ UnitC=100; % Change of meter into centimeter; Hc=0.02; % Henry's constant; GWT=7; % The gain factor(dimensionless),which assesses the temperature - % dependence of the soil water retention curve is set as 7 for + % dependence of the soil water retention curve is set as 7 for % sand (Noborio et al, 1996); -MU_a=1.846*10^(-4); % (g.cm^-1.s^-1)Viscosity of air (original 1.846*10^(-5)kg.m^-1.s^-1); +MU_a=1.846*10^(-4); % (g.cm^-1.s^-1)Viscosity of air (original 1.846*10^(-5)kg.m^-1.s^-1); Gamma0=71.89; % The surface tension of soil water at 25 Cels degree. (g.s^-2); Gamma_w=RHOL*g; % Specific weight of water(g.cm^-2.s^-2); Lambda1=0.228/UnitC;%-0.197/UnitC;% 0.243/UnitC; % Coefficients in thermal conductivity; @@ -421,32 +420,38 @@ Gsc=1360; % The solar constant (1360 W��m^-2) Sigma_E=4.90*10^(-9); % The stefan-Boltzman constant.(=4.90*10^(-9) MJ��m^-2��Cels^-4��d^-1) P_g0=95197.850;%951978.50; % The mean atmospheric pressure (Should be given in new simulation period subroutine.) -rroot=1.5*1e-3; +rroot=1.5*1e-3; RTB=1000; %initial root total biomass (g m-2) Precipp=0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Input for producing initial soil moisture and soil temperature profile %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -Mdata=textread([InputPath, 'Mdata.txt']); -Ta_msr=Mdata(:,2)'; +fID = fopen([InputPath, 'Mdata.txt']); +Mdata = textscan(fID, ""); +fclose(fID); + +Ta_msr=Mdata{:,2}; +RH_msr=Mdata{:,3}; +WS_msr=Mdata{:,4}; +Pg_msr=Mdata{:,5}; +Rn_msr=Mdata{:,7}; +Rns_msr=Mdata{:,7}; +VPD_msr=Mdata{:,9}; +LAI_msr=Mdata{:,10}; +G_msr=Mdata{:,7}*0.15; +Ts_msr=Mdata{:,2}; +Precip_msr=Mdata{:,6}*10*DELT; + for jj=1:Dur_tot -if Ta_msr(jj)<-100 -Ta_msr(jj)=NaN; -end + if Ta_msr(jj)<-100 + Ta_msr(jj)=NaN; + end end -RH_msr=Mdata(:,3)'; -WS_msr=Mdata(:,4)'; -Pg_msr=Mdata(:,5)'; -Rn_msr=Mdata(:,7)'; -Rns_msr=Mdata(:,7)'; -VPD_msr=Mdata(:,9)'; -LAI_msr=Mdata(:,10)'; -G_msr=Mdata(:,7)'*0.15; -Ts_msr=Mdata(:,2)'; -Precip_msr=Mdata(:,6)'*10*DELT; Tmin=min(Ts_msr); + Precip_msr=Precip_msr.*(1-fmax.*exp(-0.5*0.5*Tot_Depth/100)); + % load initial soil moisture and soil temperature from ERA5 Initial_path=dir(fullfile(InitialConditionPath,[sitename,'*.nc'])); InitND1=5; % Unit of it is cm. These variables are used to indicated the depth corresponding to the measurement. @@ -455,7 +460,7 @@ InitND4=100; InitND5=200; InitND6=300; -if SWCC==0 +if SWCC==0 InitT0= -1.762; %-1.75estimated soil surface temperature-1.762 InitT1= -0.662; InitT2= 0.264; @@ -464,7 +469,7 @@ InitT5= 3.657;%; InitT6= 6.033; BtmT=6.62; %9 8.1 -InitX0= 0.088; +InitX0= 0.088; InitX1= 0.095; % Measured soil liquid moisture content InitX2= 0.180; %0.169 InitX3= 0.213; %0.205 @@ -481,7 +486,7 @@ InitT5= ncread([InitialConditionPath,Initial_path(5).name],'stl4')-273.15; InitT6= ncread([InitialConditionPath,Initial_path(5).name],'stl4')-273.15; Tss = InitT0; - if InitT0 < 0 || InitT1 < 0 || InitT2 < 0 || InitT3 < 0 || InitT4 < 0 || InitT5 < 0 || InitT6 < 0 + if InitT0 < 0 || InitT1 < 0 || InitT2 < 0 || InitT3 < 0 || InitT4 < 0 || InitT5 < 0 || InitT6 < 0 InitT0= 0; InitT1= 0; InitT2= 0; @@ -513,6 +518,6 @@ InitX4= fieldMC(4); %0.14335 InitX5= fieldMC(5); InitX6= fieldMC(6); -BtmX = fieldMC(6); +BtmX = fieldMC(6); end end \ No newline at end of file diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index de3d57f2..f2e676d3 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -1,12 +1,12 @@ %% STEMMUS-SCOPE.m (script) -% STEMMUS-SCOPE is a model for Integrated modeling of canopy photosynthesis, fluorescence, -% and the transfer of energy, mass, and momentum in the soil–plant–atmosphere continuum +% STEMMUS-SCOPE is a model for Integrated modeling of canopy photosynthesis, fluorescence, +% and the transfer of energy, mass, and momentum in the soil–plant–atmosphere continuum % % Version: 1.0.1 % % Copyright (C) 2021 Yunfei Wang, Lianyu Yu, Yijian Zeng, Christiaan Van der Tol, Bob Su -% Contact: y.wang-3@utwente.nl; l.yu@utwente.nl; y.zeng@utwente.nl; c.vandertol@utwente.nl; z.su@utwente.nl +% Contact: y.wang-3@utwente.nl; l.yu@utwente.nl; y.zeng@utwente.nl; c.vandertol@utwente.nl; z.su@utwente.nl % % This program is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by @@ -26,37 +26,50 @@ % but there still global variables here, because we not sure which % progresses related to these global variables. +% Load in required Octave packages if STEMMUS-SCOPE is being run in Octave: +if exist('OCTAVE_VERSION', 'builtin') ~= 0 + pkg load netcdf + pkg load statistics +end + % Read the configPath file. Due to using MATLAB compiler, we cannot use run(CFG) global CFG if isempty(CFG) CFG = '../config_file_crib.txt'; end -disp (['Reading config from ',CFG]) -[DataPaths.soilProperty, DataPaths.input, DataPaths.output, ... - DataPaths.forcingPath, forcingFileName, numberOfTimeSteps, ... - DataPaths.initialCondition] = io.read_config(CFG); +disp (['Reading config from ', CFG]) +global InputPath OutputPath InitialConditionPath +[InputPath, OutputPath, InitialConditionPath] = io.read_config(CFG); -% Prepare forcing data +% Prepare forcing and soil data global IGBP_veg_long latitude longitude reference_height canopy_height sitename DELT Dur_tot -[SiteProperties, DELT, forcingTimeLength] = io.prepareForcingData(DataPaths, forcingFileName); -SoilPropertyPath = DataPaths.soilProperty; -InputPath = DataPaths.input; -OutputPath = DataPaths.output; -InitialConditionPath = DataPaths.initialCondition; -IGBP_veg_long = SiteProperties.igbpVegLong; +[SiteProperties, SoilProperties, TimeProperties] = io.prepareInputData(InputPath); +IGBP_veg_long = SiteProperties.IGBP_veg_long; latitude = SiteProperties.latitude; longitude = SiteProperties.longitude; -reference_height = SiteProperties.referenceHeight; -canopy_height = SiteProperties.canopyHeight; -sitename = SiteProperties.siteName; - -%Set the end time of the main loop in STEMMUS_SCOPE.m -%using config file or time length of forcing file -if isnan(numberOfTimeSteps) - Dur_tot=forcingTimeLength; -else - Dur_tot = min(numberOfTimeSteps, forcingTimeLength); -end +reference_height = SiteProperties.reference_height; +canopy_height = SiteProperties.canopy_height; +sitename = SiteProperties.sitename; + +DELT = TimeProperties.DELT; +Dur_tot = TimeProperties.Dur_tot; + +global SaturatedK SaturatedMC ResidualMC Coefficient_n Coefficient_Alpha +global porosity FOC FOS MSOC Coef_Lamda fieldMC fmax theta_s0 Ks0 +SaturatedK = SoilProperties.SaturatedK; +SaturatedMC = SoilProperties.SaturatedMC; +ResidualMC = SoilProperties.ResidualMC; +Coefficient_n = SoilProperties.Coefficient_n; +Coefficient_Alpha = SoilProperties.Coefficient_Alpha; +porosity = SoilProperties.porosity; +FOC = SoilProperties.FOC; +FOS = SoilProperties.FOS; +MSOC = SoilProperties.MSOC; +Coef_Lamda = SoilProperties.Coef_Lamda; +fieldMC = SoilProperties.fieldMC; +fmax = SoilProperties.fmax; +theta_s0 = SoilProperties.theta_s0; +Ks0 = SoilProperties.Ks0; %% run Constants %input soil parameters @@ -72,8 +85,8 @@ global SAVEKfL_h KCHK FCHK TKCHK hCHK TSAVEhh SAVEhh_frez Precip SAVETT INDICATOR thermal global Theta_g MU_a DeltZ Alpha_Lg global Beta_g KaT_Switch -global D_V D_A fc Eta nD POR -global ThmrlCondCap ZETA XK DVT_Switch +global D_V D_A fc Eta nD POR +global ThmrlCondCap ZETA XK DVT_Switch global MU_W Ks RHOL global Lambda1 Lambda2 Lambda3 RHO_bulk global RHODA RHOV c_a c_V c_L @@ -83,7 +96,7 @@ global hThmrl Tr COR IS Hystrs XWRE global Theta_V IH CORh XCAP Phi_s RHOI Lamda global W WW D_Ta SSUR -global W_Chg SWCC Imped Ratio_ice ThermCond Beta_gBAR Alpha_LgBAR +global W_Chg SWCC Imped Ratio_ice ThermCond Beta_gBAR Alpha_LgBAR global KLT_Switch xERR hERR TERR NBChB NBCT uERR global L TCON_dry TPS1 TPS2 TCON0 XSOC TCON_s global HCAP SF TCA GA1 GA2 GB1 GB2 HCD ZETA0 CON0 PS1 PS2 XWILT FEHCAP QMTT QMBB Evapo trap RnSOIL PrecipO @@ -102,12 +115,12 @@ if useXLSX == 0 % parameter_file = { 'setoptions.m', 'filenames.m', 'inputdata.txt'}; % Read parameter file which is 'input_data.xlsx' and return it as options. -% options = io.setOptions(parameter_file,path_input); +% options = io.setOptions(parameter_file,path_input); warning("the current stemmus-scope does not support useXLSX=0"); else - parameter_file = {'input_data.xlsx'}; - options = io.readStructFromExcel([path_input char(parameter_file)], 'options', 3, 1); -end + parameter_file = {'input_data.xlsx'}; + options = io.readStructFromExcel([path_input char(parameter_file)], 'options', 2, 1); +end if options.simulation>2 || options.simulation<0, fprintf('\n simulation option should be between 0 and 2 \r'); return, end @@ -116,7 +129,7 @@ if useXLSX==0 run([path_input parameter_file{2}]) else - [dummy,X] = xlsread([path_input char(parameter_file)],'filenames'); + [dummy,X] = xlsread([path_input char(parameter_file)],'filenames'); j = find(~strcmp(X(:,2),{''})); X = X(j,(1:end)); end @@ -154,7 +167,7 @@ 'I will use 0.25*Cab instead'); options.Cca_function_of_Cab = 1; else - + if ~(options.simulation==1) && (i==30 || i==32) warning('warning: input "', V(i).Name, '" not provided in input spreadsheet...', ... 'I will use the MODTRAN spectrum as it is'); @@ -185,7 +198,7 @@ end end end - + % the current stemmus-scope does not support useXLSX=0 if useXLSX==0 j2 = []; j1 = j+1; @@ -199,8 +212,8 @@ else V(i).Val = N(j2); end - - + + else if sum(~isnan(N(j,:)))<1 V(i).Val = -999; @@ -218,64 +231,64 @@ V(55).Val=mean(Ta_msr); %Input T parameters for different vegetation type sitename1=cellstr(sitename); -if strcmp(IGBP_veg_long(1:18)', 'Permanent Wetlands') +if strcmp(IGBP_veg_long(1:18)', 'Permanent Wetlands') V(14).Val= [0.2 0.3 288 313 328]; % These are five parameters specifying the temperature response. V(9).Val= [120]; % Vcmax, maximum carboxylation capacity (at optimum temperature) V(10).Val= [9]; % Ball-Berry stomatal conductance parameter V(11).Val= [0]; % Photochemical pathway: 0=C3, 1=C4 V(28).Val= [0.05]; % leaf width -elseif strcmp(IGBP_veg_long(1:19)', 'Evergreen Broadleaf') +elseif strcmp(IGBP_veg_long(1:19)', 'Evergreen Broadleaf') V(14).Val= [0.2 0.3 283 311 328]; V(9).Val= [80]; V(10).Val= [9]; V(11).Val= [0]; V(28).Val= [0.05]; -elseif strcmp(IGBP_veg_long(1:19)', 'Deciduous Broadleaf') +elseif strcmp(IGBP_veg_long(1:19)', 'Deciduous Broadleaf') V(14).Val= [0.2 0.3 283 311 328]; V(9).Val= [80]; V(10).Val= [9]; - V(11).Val= [0]; + V(11).Val= [0]; V(28).Val= [0.05]; -elseif strcmp(IGBP_veg_long(1:13)', 'Mixed Forests') +elseif strcmp(IGBP_veg_long(1:13)', 'Mixed Forests') V(14).Val= [0.2 0.3 281 307 328]; V(9).Val= [80]; V(10).Val= [9]; - V(11).Val= [0]; + V(11).Val= [0]; V(28).Val= [0.04]; -elseif strcmp(IGBP_veg_long(1:20)', 'Evergreen Needleleaf') +elseif strcmp(IGBP_veg_long(1:20)', 'Evergreen Needleleaf') V(14).Val= [0.2 0.3 278 303 328]; V(9).Val= [80]; V(10).Val= [9]; - V(11).Val= [0]; + V(11).Val= [0]; V(28).Val= [0.01]; -elseif strcmp(IGBP_veg_long(1:9)', 'Croplands') +elseif strcmp(IGBP_veg_long(1:9)', 'Croplands') if isequal(sitename1,{'ES-ES2'})||isequal(sitename1,{'FR-Gri'})||isequal(sitename1,{'US-ARM'})||isequal(sitename1,{'US-Ne1'}) V(14).Val= [0.2 0.3 278 303 328]; V(9).Val= [50]; V(10).Val= [4]; - V(11).Val= [1]; + V(11).Val= [1]; V(13).Val= [0.025]; % Respiration = Rdparam*Vcmcax V(28).Val= [0.03]; - else + else V(14).Val= [0.2 0.3 278 303 328]; V(9).Val= [120]; V(10).Val= [9]; - V(11).Val= [0]; - V(28).Val= [0.03]; + V(11).Val= [0]; + V(28).Val= [0.03]; end elseif strcmp(IGBP_veg_long(1:15)', 'Open Shrublands') V(14).Val= [0.2 0.3 288 313 328]; V(9).Val= [120]; V(10).Val= [9]; - V(11).Val= [0]; + V(11).Val= [0]; V(28).Val= [0.05]; -elseif strcmp(IGBP_veg_long(1:17)', 'Closed Shrublands') +elseif strcmp(IGBP_veg_long(1:17)', 'Closed Shrublands') V(14).Val= [0.2 0.3 288 313 328]; V(9).Val= [80]; V(10).Val= [9]; V(11).Val= [0]; V(28).Val= [0.05]; -elseif strcmp(IGBP_veg_long(1:8)', 'Savannas') +elseif strcmp(IGBP_veg_long(1:8)', 'Savannas') V(14).Val= [0.2 0.3 278 313 328]; V(9).Val= [120]; V(10).Val= [9]; @@ -287,13 +300,13 @@ V(10).Val= [9]; V(11).Val= [0]; V(28).Val= [0.03]; -elseif strcmp(IGBP_veg_long(1:9)', 'Grassland') +elseif strcmp(IGBP_veg_long(1:9)', 'Grassland') V(14).Val= [0.2 0.3 288 303 328]; if isequal(sitename1,{'AR-SLu'})||isequal(sitename1,{'AU-Ync'})||isequal(sitename1,{'CH-Oe1'})||isequal(sitename1,{'DK-Lva'})||isequal(sitename1,{'US-AR1'})||isequal(sitename1,{'US-AR2'})||isequal(sitename1,{'US-Aud'})||isequal(sitename1,{'US-SRG'}) V(9).Val= [120]; V(10).Val= [4]; V(11).Val= [1]; - V(13).Val= [0.025]; + V(13).Val= [0.025]; else V(9).Val= [120]; V(10).Val= [9]; @@ -417,16 +430,16 @@ atmo.M = helpers.aggreg(atmfile,spectral.SCOPEspec); %% 13. create output files -[Output_dir, f, fnames] = io.create_output_files_binary(parameter_file, sitename, path_of_code, path_input, path_output, spectral, options); +[Output_dir, fnames] = io.create_output_files_binary(parameter_file, sitename, path_of_code, path_input, path_output, spectral, options); run StartInit; % Initialize Temperature, Matric potential and soil air pressure. %% 14. Run the model -fprintf('\n The calculations start now \r') +disp('The calculations start now') calculate = 1; TIMEOLD=0;SAVEhh_frez=zeros(NN+1,1);FCHK=zeros(1,NN);KCHK=zeros(1,NN);hCHK=zeros(1,NN); TIMELAST=0; -SAVEtS=tS; kk=0; %DELT=Delt_t; +SAVEtS=tS; kk=0; %DELT=Delt_t; for i = 1:1:Dur_tot KT=KT+1; % Counting Number of timesteps if KT>1 && Delt_t>(TEND-TIME) @@ -454,7 +467,7 @@ TOLD_CRIT(MN)=T_CRIT(MN); T_CRIT(MN)=TT_CRIT(MN); %TTT_CRIT(MN,KT)=TT_CRIT(MN); - + hOLD(MN)=h(MN); h(MN)=hh(MN); %hhh(MN,KT)=hh(MN); @@ -503,20 +516,20 @@ else calculate = 1; end - + if calculate - + iter.counter = 0; - + LIDF_file = char(F(22).FileName); if ~isempty(LIDF_file) canopy.lidf = dlmread([path_input,'leafangles/',LIDF_file],'',3,0); else canopy.lidf = equations.leafangles(canopy.LIDFa,canopy.LIDFb); % This is 'ladgen' in the original SAIL model, end - - leafbio.emis = 1-leafbio.rho_thermal-leafbio.tau_thermal; - + + leafbio.emis = 1-leafbio.rho_thermal-leafbio.tau_thermal; + if options.calc_PSI fversion = @fluspect_B_CX; else @@ -526,15 +539,15 @@ leafopt = fversion(spectral,leafbio,optipar); leafbio.V2Z = 1; leafoptZ = fversion(spectral,leafbio,optipar); - + IwlP = spectral.IwlP; IwlT = spectral.IwlT; - + rho(IwlP) = leafopt.refl; tau(IwlP) = leafopt.tran; rlast = rho(nwlP); tlast = tau(nwlP); - + if options.soilspectrum == 0 rs(IwlP) = rsfile(:,soil.spectrum+1); else @@ -543,7 +556,7 @@ rs(IwlP) = BSM(soil,optipar,soilemp); end rslast = rs(nwlP); - + switch options.rt_thermal case 0 rho(IwlT) = ones(nwlT,1) * leafbio.rho_thermal; @@ -556,26 +569,26 @@ end leafopt.refl = rho; % extended wavelength ranges are stored in structures leafopt.tran = tau; - + reflZ = leafopt.refl; tranZ = leafopt.tran; reflZ(1:300) = leafoptZ.refl(1:300); tranZ(1:300) = leafoptZ.tran(1:300); leafopt.reflZ = reflZ; leafopt.tranZ = tranZ; - + soil.refl = rs; - + soil.Ts = meteo.Ta * ones(2,1); % initial soil surface temperature - + if length(F(4).FileName)>1 && options.simulation==0 atmfile = [path_input 'radiationdata/' char(F(4).FileName(k))]; atmo.M = helpers.aggreg(atmfile,spectral.SCOPEspec); end atmo.Ta = meteo.Ta; - + [rad,gap,profiles] = RTMo(spectral,atmo,soil,leafopt,canopy,angles,meteo,rad,options); - + switch options.calc_ebal case 1 [iter,fluxes,rad,thermal,profiles,soil,RWU,frac] ... @@ -591,15 +604,15 @@ if options.calc_xanthophyllabs [rad] = RTMz(spectral,rad,soil,leafopt,canopy,gap,angles,profiles); end - + if options.calc_planck rad = RTMt_planck(spectral,rad,soil,leafopt,canopy,gap,angles,thermal.Tcu,thermal.Tch,thermal.Ts(2),thermal.Ts(1),1); end - + if options.calc_directional directional = calc_brdf(options,directional,spectral,angles,rad,atmo,soil,leafopt,canopy,meteo,profiles,thermal); end - + otherwise Fc = (1-gap.Ps(1:end-1))'/nl; % Matrix containing values for Ps of canopy fluxes.aPAR = canopy.LAI*(Fc*rad.Pnh + equations.meanleaf(canopy,rad.Pnu , 'angles_and_layers',gap.Ps));% net PAR leaves @@ -621,42 +634,42 @@ else rad.Femtot = 1E3*leafbio.fqe* optipar.phi(spectral.IwlF) * fluxes.aPAR_Cab_eta; end - end + end end if options.simulation==2 && telmax>1, vi = helpers.count(nvars,vi,vmax,1); end - if KT==1 + if KT==1 if isreal(fluxes.Actot)&&isreal(thermal.Tsave)&&isreal(fluxes.lEstot)&&isreal(fluxes.lEctot) Acc=fluxes.Actot; lEstot =fluxes.lEstot; lEctot =fluxes.lEctot; - Tss=thermal.Tsave; + Tss=thermal.Tsave; else Acc=0; lEstot =0; lEctot =0; - Tss=Ta_msr(KT); + Tss=Ta_msr(KT); end elseif NoTime(KT)>NoTime(KT-1) if isreal(fluxes.Actot)&&isreal(thermal.Tsave)&&isreal(fluxes.lEstot)&&isreal(fluxes.lEctot) Acc=fluxes.Actot; lEstot =fluxes.lEstot; lEctot =fluxes.lEctot; - Tss=thermal.Tsave; + Tss=thermal.Tsave; else Acc=0; lEstot =0; lEctot =0; - Tss=Ta_msr(KT); + Tss=Ta_msr(KT); end end - + sfactortot(KT)=sfactor; - DSTOR0=DSTOR; - + DSTOR0=DSTOR; + if KT>1 run SOIL1 end - + for ML=1:NL %QVV(ML,KT)=QV(ML); %QLL(ML,KT)=QL(ML,1); @@ -729,7 +742,7 @@ end end run Forcing_PARM - + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for KIT=1:NIT % Start the iteration procedure in a time step. [TT_CRIT,hh_frez]=HT_frez(hh,T0,g,L_f,TT,NN,hd,Tmin); @@ -760,15 +773,15 @@ DSTOR=min(EXCESS,DSTMAX); RS(KT)=(EXCESS-DSTOR)/Delt_t; end - + if Soilairefc==1 run Air_sub; end - + if Thmrlefc==1 run Enrgy_sub; end - + if max(CHK)<0.1 %&& max(FCHK)<0.001 %&& max(hCHK)<0.001 %&& min(KCHK)>0.001 break end @@ -806,7 +819,7 @@ hOLD(MN)=h(MN); h(MN)=hh(MN); %hhh(MN,KT)=hh(MN); - %HRA(MN,KT)=HR(MN); + %HRA(MN,KT)=HR(MN); if Thmrlefc==1 TOLD(MN)=T(MN); T(MN)=TT(MN); @@ -838,9 +851,14 @@ %SAVEDSTOR(KT)=DSTOR; end kk=k; - n_col = io.output_data_binary(f, k, xyt, rad, canopy, V, vi, vmax, options, fluxes, meteo, iter, thermal, spectral, gap, profiles, Sim_Theta_U, Sim_Temp, Trap, Evap); + + % Open files for writing + file_ids = structfun(@(x) fopen(x, 'a'), fnames, 'UniformOutput',false); + n_col = io.output_data_binary(file_ids, k, xyt, rad, canopy, V, vi, vmax, options, fluxes, meteo, iter, thermal, spectral, gap, profiles, Sim_Theta_U, Sim_Temp, Trap, Evap); + fclose("all"); end -fprintf('\n The calculations end now \r') + +disp('The calculations end now') if options.verify io.output_verification(Output_dir) end @@ -855,5 +873,5 @@ save([Output_dir,'output.mat']) %if options.makeplots % plot.plots(Output_dir) -%end +%end diff --git a/src/debug_Octave.m b/src/debug_Octave.m new file mode 100644 index 00000000..2186ee36 --- /dev/null +++ b/src/debug_Octave.m @@ -0,0 +1,2 @@ +config_file = "C:/STEMMUS_SCOPE_data/workdir/input/ZA-Kru_2022-09-30-1448/ZA-Kru_2022-09-30-1448_config.txt" +STEMMUS_SCOPE_exe(config_file); diff --git a/src/ebal.m b/src/ebal.m index d9a4289d..3cee7b8c 100644 --- a/src/ebal.m +++ b/src/ebal.m @@ -1,4 +1,4 @@ -function [iter,fluxes,rad,thermal,profiles,soil,RWU,frac] ... +function [iter,fluxes,rad,thermal,profiles,soil,RWU,frac] ... = ebal(iter,options,spectral,rad,gap,leafopt, ... angles,meteo,soil,canopy,leafbio,xyt,k,profiles,Delt_t) global Rl DeltZ Ks Theta_s Theta_r Theta_LL bbx NL KT sfactor PSItot sfactortot Theta_f @@ -27,7 +27,7 @@ % soil_respiration.m % % Table of contents of the function -% +% % 1. Initialisations for the iteration loop % intial values are attributed to variables % 2. Energy balance iteration loop @@ -38,29 +38,29 @@ % % The energy balance iteration loop works as follows: % -% RTMo More or less the classic SAIL model for Radiative +% RTMo More or less the classic SAIL model for Radiative % Transfer of sun and sky light (no emission by the vegetation) -% While continue Here an iteration loop starts to close the energy -% balance, i.e. to match the micro-meteorological model +% While continue Here an iteration loop starts to close the energy +% balance, i.e. to match the micro-meteorological model % and the radiative transfer model -% RTMt_sb A numerical Radiative Transfer Model for thermal +% RTMt_sb A numerical Radiative Transfer Model for thermal % radiation emitted by the vegetation -% resistances Calculates aerodynamic and boundary layer resistances +% resistances Calculates aerodynamic and boundary layer resistances % of vegetation and soil (the micro-meteorological model) -% biochemical Calculates photosynthesis, fluorescence and stomatal +% biochemical Calculates photosynthesis, fluorescence and stomatal % resistance of leaves (or biochemical_MD12: alternative) -% heatfluxes Calculates sensible and latent heat flux of soil and +% heatfluxes Calculates sensible and latent heat flux of soil and % vegetation -% Next soil heat flux is calculated, the energy balance -% is evaluated, and soil and leaf temperatures adjusted +% Next soil heat flux is calculated, the energy balance +% is evaluated, and soil and leaf temperatures adjusted % to force energy balance closure % end {while continue} -% -% meanleaf Integrates the fluxes over all leaf inclinations +% +% meanleaf Integrates the fluxes over all leaf inclinations % azimuth angles and layers, integrates over the spectrum % % usage: -%[iter,fluxes,rad,profiles,thermal] ... +%[iter,fluxes,rad,profiles,thermal] ... % = ebal(iter,options,spectral,rad,gap,leafopt, ... % angles,meteo,soil,canopy,leafbio) % @@ -68,7 +68,7 @@ % specified in a readme file. % % Input: -% +% % iter numerical parameters used in the iteration for energy balance closure % options calculation options % spectral spectral resolutions and wavelengths @@ -81,7 +81,7 @@ % leafbio leaf biochemical parameters % % Output: -% +% % iter numerical parameters used in the iteration for energy balance closure % fluxes energy balance, turbulent, and CO2 fluxes % rad radiation spectra @@ -100,7 +100,7 @@ CONT = 1; % is 0 when the calculation has finished es_fun = @(T)6.107*10.^(7.5.*T./(237.3+T)); s_fun = @(es, T) es*2.3026*7.5*237.3./(237.3+T).^2; -t = xyt.t(k); +t = xyt.t(k); Ta = meteo.Ta; ea = meteo.ea; Ca = meteo.Ca; @@ -112,7 +112,7 @@ else Deltat = Delt_t; end - x = [1:12;1:12]'*Deltat; + x = [1:12;1:12]'.*Deltat; Tsold = soil.Tsold; end @@ -174,31 +174,31 @@ %'Energy balance loop (Energy balance and radiative transfer) while CONT % while energy balance does not close - + % 2.1. Net radiation of the components % Thermal radiative transfer model for vegetation emission (with Stefan-Boltzman's equation) rad = RTMt_sb(spectral,rad,soil,leafopt,canopy,gap,angles,Tcu,Tch,Ts(2),Ts(1),1); % Add net radiation of (1) solar and sky and (2) thermal emission model - + Rnhct = rad.Rnhct; Rnuct = rad.Rnuct; Rnhst = rad.Rnhst; Rnust = rad.Rnust; - + Rnhc = rad.Rnhc; Rnuc = rad.Rnuc; Rnhs = rad.Rnhs; Rnus = rad.Rnus; - + Rnch = Rnhc + Rnhct; % Canopy (shaded) net radiation Rncu = Rnuc + Rnuct; % Canopy (sunlit) net radiation Rnsh = Rnhs + Rnhst; % Soil (shaded) net radiation Rnsu = Rnus + Rnust; % Soil (sunlit) net radiation Rns = [Rnsh Rnsu]'; % Soil (sun+sh) net radiation - + % 2.2. Aerodynamic roughness % calculate friction velocity [m s-1] and aerodynamic resistances [s m-1] - + resist_in.u = max(meteo.u,.2); resist_in.L = L; resist_in.LAI = canopy.LAI; @@ -211,16 +211,16 @@ resist_in.hc = canopy.hc; resist_in.w = canopy.leafwidth; resist_in.Cd = canopy.Cd; - + [resist_out] = resistances(resist_in); - + ustar = resist_out.ustar; raa = resist_out.raa; rawc = resist_out.rawc; raws = resist_out.raws; - + % 2.3. Biochemical processes - + % photosynthesis (A), fluorescence factor (F), and stomatal resistance (rcw), for shaded (1) and sunlit (h) leaves biochem_in.Fluorescence_model = options.Fluorescence_model; biochem_in.Type = leafbio.Type; @@ -229,7 +229,7 @@ biochem_in.BallBerry0 = leafbio.BallBerry0; biochem_in.O = meteo.Oa; biochem_in.Rdparam = leafbio.Rdparam; - + if options.Fluorescence_model==2 % specific for the v.Caemmerer-Magnani model b = @biochemical_MD12; biochem_in.Tyear = leafbio.Tyear; @@ -241,7 +241,7 @@ b = @biochemical; % specific for Berry-v.d.Tol model biochem_in.tempcor = options.apply_T_corr; biochem_in.Tparams = leafbio.Tparam; - biochem_in.stressfactor = SMCsf; + biochem_in.stressfactor = SMCsf; end % for shaded leaves @@ -250,7 +250,7 @@ biochem_in.Vcmo = fVh.*leafbio.Vcmo; biochem_in.Cs = Cch; biochem_in.Q = rad.Pnh_Cab*1E6; - + biochem_out = b(biochem_in); Ah = biochem_out.A; Ahh = biochem_out.Ag; @@ -259,16 +259,16 @@ rcwh = biochem_out.rcw; qEh = biochem_out.qE; % vCaemmerer- Magnani does not generate this parameter (dummy value) Knh = biochem_out.Kn; - + % for sunlit leaves biochem_in.T = Tcu; biochem_in.eb = ecu; biochem_in.Vcmo = fVu.*leafbio.Vcmo; biochem_in.Cs = Ccu; biochem_in.Q = rad.Pnu_Cab*1E6; - + biochem_out = b(biochem_in); - + Au = biochem_out.A; % Ag? or A? Auu = biochem_out.Ag; %GPP calculation. Ciu = biochem_out.Ci; @@ -276,17 +276,17 @@ rcwu = biochem_out.rcw; qEu = biochem_out.qE; Knu = biochem_out.Kn; - + Pinh = rad.Pnh; Pinu = rad.Pnu; Pinh_Cab = rad.Pnh_Cab; Pinu_Cab = rad.Pnu_Cab; - Rnh_PAR = rad.Rnh_PAR; + Rnh_PAR = rad.Rnh_PAR; Rnu_PAR = rad.Rnu_PAR; - + % 2.4. Fluxes (latent heat flux (lE), sensible heat flux (H) and soil heat flux G % in analogy to Ohm's law, for canopy (c) and soil (s). All in units of [W m-2] - + %soil.PSIs; rss = soil.rss; rac = (LAI+1)*(raa+rawc); @@ -295,7 +295,7 @@ [lEch,Hch,ech,Cch,lambdah,sh] = heatfluxes(rac,rcwh,Tch,ea,Ta,e_to_q,PSI,Ca,Cih,constants,es_fun,s_fun); [lEcu,Hcu,ecu,Ccu,lambdau,su] = heatfluxes(rac,rcwu,Tcu,ea,Ta,e_to_q,PSI,Ca,Ciu,constants,es_fun,s_fun); [lEs,Hs,~,~,lambdas,ss] = heatfluxes(ras,rss,Ts ,ea,Ta,e_to_q,PSIss,Ca,Ca,constants,es_fun,s_fun); - + %if any( ~isreal( Cch )) || any( ~isreal( Ccu(:) )) % error('Heatfluxes produced complex values for CO2 concentration!') %end @@ -322,7 +322,7 @@ BB2=AA2(~isinf(AA2)); PSI1 = (sum(BB1)-Trans)/sum(BB2); if isnan(PSI1) - PSI1 = -1; + PSI1 = -1; end if ~isreal(PSI1) PSI1 = -1; @@ -336,49 +336,49 @@ %%%%%%% if SoilHeatMethod==2 G = 0.30*Rns; - else - G = GAM/sqrt(pi)*2*sum(([Ts'; Tsold(1:end-1,:)] - Tsold)/Deltat .* (sqrt(x) - sqrt(x-Deltat))); + else + G = GAM/sqrt(pi)*2*sum(([Ts'; Tsold(1:end-1,:)] - Tsold)./Deltat .* (sqrt(x) - sqrt(x-Deltat))); G = G'; end % 2.5. Monin-Obukhov length L L = -rhoa*cp*ustar.^3.*(Ta+273.15)./(kappa*g*Htot); % [1] - L(L<-1E3) = -1E3; % [1] - L(L>1E2) = 1E2; % [1] - L(isnan(L)) = -1; % [1] - + L(L<-1E3) = -1E3; % [1] + L(L>1E2) = 1E2; % [1] + L(isnan(L)) = -1; % [1] + % 2.6. energy balance errors, continue criterion and iteration counter EBerch = Rnch -lEch -Hch; EBercu = Rncu -lEcu -Hcu; EBers = Rns -lEs -Hs - G; - + counter = counter+1; % Number of iterations maxEBercu = max(max(max(abs(EBercu)))); maxEBerch = max(abs(EBerch)); maxEBers = max(abs(EBers)); - - CONT = ( maxEBercu > maxEBer |... + + CONT = ( maxEBercu > maxEBer |... maxEBerch > maxEBer |... maxEBers > maxEBer) &... counter < maxit+1;% Continue iteration? if counter==5, Wc = 0.6; end if counter==10; Wc = 0.4; end - if counter==30; Wc = 0.1; end - - % 2.7. New estimates of soil (s) and leaf (c) temperatures, shaded (h) and sunlit (1) + if counter==30; Wc = 0.1; end + + % 2.7. New estimates of soil (s) and leaf (c) temperatures, shaded (h) and sunlit (1) Tch = Tch + Wc*EBerch./((rhoa*cp)./rac + rhoa*lambdah*e_to_q.*sh./(rac+rcwh)+ 4*leafbio.emis*sigmaSB*(Tch+273.15).^3); Tcu = Tcu + Wc*EBercu./((rhoa*cp)./rac + rhoa*lambdau*e_to_q.*su./(rac+rcwu)+ 4*leafbio.emis*sigmaSB*(Tcu+273.15).^3); - Ts = Ts + Wc*EBers./(rhoa*cp./ras + rhoa*lambdas*e_to_q.*ss/(ras+rss)+ 4*(1-soil.rs_thermal)*sigmaSB*(Ts+273.15).^3); % Ts contains shaded soil temperature and sunlit soil temperature + Ts = Ts + Wc*EBers./(rhoa*cp./ras + rhoa*lambdas*e_to_q.*ss/(ras+rss)+ 4*(1-soil.rs_thermal)*sigmaSB*(Ts+273.15).^3); % Ts contains shaded soil temperature and sunlit soil temperature Tch(abs(Tch)>100) = Ta; Tcu(abs(Tcu)>100) = Ta; Ts(abs(Ts)>100) = Ta; if (any(isnan(Tch)) || any(isnan(Tcu(:)))), warning('Canopy temperature gives NaNs'), end if any(isnan(Ts)), warning('Soil temperature gives NaNs'), end - + end iter.counter = counter; profiles.etah = Fh; -profiles.etau = Fu; +profiles.etau = Fu; if SoilHeatMethod<2 Tsold(2:end,:) = soil.Tsold(1:end-1,:); @@ -402,11 +402,11 @@ '\n Energy balance error sunlit vegetation = %4.2f W m-2 ', ... '\n Energy balance error shaded vegetation = %4.2f W m-2 ', ... '\n Energy balance error soil = %4.2f W m-2 '],maxEBercu, maxEBerch, maxEBers); - + end %% 4. Calculate the output per layer -if options.calc_vert_profiles +if options.calc_vert_profiles [Hcu1d ] = equations.meanleaf(canopy,Hcu, 'angles'); % [nli,nlo,nl] mean sens heat sunlit leaves [lEcu1d ] = equations.meanleaf(canopy,lEcu, 'angles'); % [nli,nlo,nl] mean latent sunlit leaves [Au1d ] = equations.meanleaf(canopy,Au, 'angles'); % [nli,nlo,nl] mean phots sunlit leaves @@ -416,7 +416,7 @@ %[Pnu1d_Cab ] = equations.meanleaf(canopy,Pinu_Cab, 'angles'); % [nli,nlo,nl] mean net radiation sunlit leaves [Rnu1d ] = equations.meanleaf(canopy,Rncu, 'angles'); % [nli,nlo,nl] mean net PAR sunlit leaves [Tcu1d ] = equations.meanleaf(canopy,Tcu, 'angles'); % [nli,nlo,nl] mean temp sunlit leaves - + profiles.Tchave = mean(Tch); % [1] mean temp shaded leaves profiles.Tch = Tch; % [nl] profiles.Tcu1d = Tcu1d; % [nl] @@ -430,10 +430,10 @@ %profiles.Pn1d_Cab = ((1-Ps(1:nl)).*Pinh_Cab + Ps(1:nl).*(Pnu1d_Cab)); %[nl] mean photos leaves, per layer profiles.Rn1d = ((1-Ps(1:nl)).*Rnch + Ps(1:nl).*(Rnu1d)); %[nl] end - + %% 5. Calculate spectrally integrated energy, water and CO2 fluxes -% sum of all leaves, and average leaf temperature +% sum of all leaves, and average leaf temperature % (note that averaging temperature is physically not correct...) Rnctot = LAI*(Fc*Rnch + equations.meanleaf(canopy,Rncu,'angles_and_layers',Ps)); % net radiation leaves @@ -445,7 +445,7 @@ Pntot = LAI*(Fc*Pinh + equations.meanleaf(canopy,Pinu,'angles_and_layers',Ps)); % net PAR leaves Pntot_Cab = LAI*(Fc*Pinh_Cab + equations.meanleaf(canopy,Pinu_Cab,'angles_and_layers',Ps)); % net PAR leaves Rntot_PAR = LAI*(Fc*Rnh_PAR + equations.meanleaf(canopy,Rnu_PAR, 'angles_and_layers',Ps));% net PAR leaves -aPAR_Cab_eta = LAI*(Fc*(profiles.etah .* Rnh_PAR) + equations.meanleaf(canopy,(profiles.etau .* Rnu_PAR), 'angles_and_layers',Ps)); +aPAR_Cab_eta = LAI*(Fc*(profiles.etah .* Rnh_PAR) + equations.meanleaf(canopy,(profiles.etau .* Rnu_PAR), 'angles_and_layers',Ps)); % ... green ePAR * relative fluorescence emission efficiency %%%%%%%%%%%%%%%%%%% [Delta_Rltot] = Root_properties(Actot,rroot); %%%%%%%%%%%%%%%%%%% Delta_Rl = fc*Delta_Rltot; @@ -506,7 +506,7 @@ thermal.Tsave = Tsave; % [oC] weighted average soil temperature thermal.raa = raa; % [s m-1] total aerodynamic resistance above canopy thermal.rawc = rawc; % [s m-1] aerodynamic resistance below canopy for canopy -thermal.raws = raws; % [s m-1] aerodynamic resistance below canopy for soil +thermal.raws = raws; % [s m-1] aerodynamic resistance below canopy for soil thermal.ustar = ustar; % [m s-1] friction velocity thermal.Tcu = Tcu; thermal.Tch = Tch; diff --git a/src/soilpropertyread.m b/src/soilpropertyread.m deleted file mode 100644 index 07e3163b..00000000 --- a/src/soilpropertyread.m +++ /dev/null @@ -1,143 +0,0 @@ -global SaturatedK SaturatedMC ResidualMC Coefficient_n Coefficient_Alpha porosity FOC FOS FOSL MSOC Coef_Lamda fieldMC latitude longitude fmax theta_s0 Ks0 sitename -% the path SoilPropertyPath is set in filereads.m -dirOutput=dir([SoilPropertyPath, 'Hydraul_Param_SoilGrids_Schaap_sl7.nc']); -%ncdisp([SoilPropertyPath,'Hydraul_Param_SoilGrids_Schaap_sl1.nc'],'/','full'); -%ncdisp([SoilPropertyPath,'Hydraul_Param_SoilGrids_Schaap_sl2.nc'],'/','full'); -%ncdisp([SoilPropertyPath,'Hydraul_Param_SoilGrids_Schaap_sl3.nc'],'/','full'); -%ncdisp([SoilPropertyPath,'Hydraul_Param_SoilGrids_Schaap_sl4.nc'],'/','full'); -%ncdisp([SoilPropertyPath,'Hydraul_Param_SoilGrids_Schaap_sl5.nc'],'/','full'); -%ncdisp([SoilPropertyPath,'Hydraul_Param_SoilGrids_Schaap_sl6.nc'],'/','full'); -%ncdisp([SoilPropertyPath,'Hydraul_Param_SoilGrids_Schaap_sl7.nc'],'/','full'); -%ncdisp([SoilPropertyPath,'CLAY1.nc','/','full'); -%% load soil property -if sitename(1:2)==['ID'] % Soil data missing at ID-Pag site, we use anothor location information here. - latitude=-1; - longitude=112; -end -lat=ncread([SoilPropertyPath,'CLAY1.nc'],'lat'); -lon=ncread([SoilPropertyPath,'CLAY1.nc'],'lon'); -for i=1:16800 - if abs(lat(i)-latitude)<0.0085 - break - end -end - for j=1:43200 - if abs(lon(j)-longitude)<0.0085 - break - end - end -depth1=ncread([SoilPropertyPath,'CLAY1.nc'],'depth'); -depth2=ncread([SoilPropertyPath,'CLAY2.nc'],'depth'); -depth3=ncread([SoilPropertyPath,'POR.nc'],'depth'); -CLAY1=ncread([SoilPropertyPath,'CLAY1.nc'],'CLAY',[j,i,1],[1,1,4]); -CLAY2=ncread([SoilPropertyPath,'CLAY2.nc'],'CLAY',[j,i,1],[1,1,4]); -SAND1=ncread([SoilPropertyPath,'SAND1.nc'],'SAND',[j,i,1],[1,1,4]); -SAND2=ncread([SoilPropertyPath,'SAND2.nc'],'SAND',[j,i,1],[1,1,4]); -SILT1=ncread([SoilPropertyPath,'SILT1.nc'],'SILT',[j,i,1],[1,1,4]); -SILT2=ncread([SoilPropertyPath,'SILT2.nc'],'SILT',[j,i,1],[1,1,4]); -OC1=ncread([SoilPropertyPath,'OC1.nc'],'OC',[j,i,1],[1,1,4]); -OC2=ncread([SoilPropertyPath,'OC2.nc'],'OC',[j,i,1],[1,1,4]); - -FOC=[CLAY1(1) CLAY1(3) CLAY2(1) CLAY2(2) CLAY2(3) CLAY2(4)]/100; %fraction of clay -FOS=[SAND1(1) SAND1(3) SAND2(1) SAND2(2) SAND2(3) SAND2(4)]/100; %fraction of sand -%FOSL=1-FOC-FOS; %fraction of silt -MSOC=double([OC1(1) OC1(3) OC2(1) OC2(2) OC2(3) OC2(4)])./10000; %mass fraction of soil organic matter -%% load lamda -lati=ncread([SoilPropertyPath,'lambda/lambda_l1.nc'],'lat'); -long=ncread([SoilPropertyPath,'lambda/lambda_l1.nc'],'lon'); -for i=1:21600 - if abs(lati(i)-latitude)<=0.0042 - break - end -end - for j=1:43200 - if abs(long(j)-longitude)<=0.0042 - break - end - end -lambda1=ncread([SoilPropertyPath,'lambda/lambda_l1.nc'],'lambda',[j,i],[1,1]); -lambda2=ncread([SoilPropertyPath,'lambda/lambda_l2.nc'],'lambda',[j,i],[1,1]); -lambda3=ncread([SoilPropertyPath,'lambda/lambda_l3.nc'],'lambda',[j,i],[1,1]); -lambda4=ncread([SoilPropertyPath,'lambda/lambda_l4.nc'],'lambda',[j,i],[1,1]); -lambda5=ncread([SoilPropertyPath,'lambda/lambda_l5.nc'],'lambda',[j,i],[1,1]); -lambda6=ncread([SoilPropertyPath,'lambda/lambda_l6.nc'],'lambda',[j,i],[1,1]); -lambda7=ncread([SoilPropertyPath,'lambda/lambda_l7.nc'],'lambda',[j,i],[1,1]); -lambda8=ncread([SoilPropertyPath,'lambda/lambda_l8.nc'],'lambda',[j,i],[1,1]); -Coef_Lamda=[lambda1 lambda3 lambda5 lambda6 lambda7 lambda8]; -%% load soil hydrulic parameters -lat=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl1_alpha.nc'],'latitude'); -lon=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl1_alpha.nc'],'longitude'); -% read data -for i=1:17924 - if abs(lat(i)-latitude)<=0.0042 - break - end -end - for j=1:43200 - if abs(lon(j)-longitude)<=0.0042 - break - end -end -% 0cm -alpha0=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl1_alpha.nc'],'alpha_0cm',[j,i],[1,1]); -n0=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl1_n.nc'],'n_0cm',[j,i],[1,1]); -theta_s0=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl1_thetas.nc'],'thetas_0cm',[j,i],[1,1]); -theta_r0=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl1_thetar.nc'],'thetar_0cm',[j,i],[1,1]); -Ks0=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl1_Ks.nc'],'Ks_0cm',[j,i],[1,1]); -% 5cm -alpha5=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl2_alpha.nc'],'alpha_5cm',[j,i],[1,1]); -n5=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl2_n.nc'],'n_5cm',[j,i],[1,1]); -theta_s5=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl2_thetas.nc'],'thetas_5cm',[j,i],[1,1]); -theta_r5=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl2_thetar.nc'],'thetar_5cm',[j,i],[1,1]); -Ks5=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl2_Ks.nc'],'Ks_5cm',[j,i],[1,1]); -% 15cm -alpha15=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl3_alpha.nc'],'alpha_15cm',[j,i],[1,1]); -n15=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl3_n.nc'],'n_15cm',[j,i],[1,1]); -theta_s15=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl3_thetas.nc'],'thetas_15cm',[j,i],[1,1]); -theta_r15=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl3_thetar.nc'],'thetar_15cm',[j,i],[1,1]); -Ks15=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl3_Ks.nc'],'Ks_15cm',[j,i],[1,1]); -% 30cm -alpha30=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl4_alpha.nc'],'alpha_30cm',[j,i],[1,1]); -n30=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl4_n.nc'],'n_30cm',[j,i],[1,1]); -theta_s30=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl4_thetas.nc'],'thetas_30cm',[j,i],[1,1]); -theta_r30=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl4_thetar.nc'],'thetar_30cm',[j,i],[1,1]); -Ks30=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl4_Ks.nc'],'Ks_30cm',[j,i],[1,1]); -% 60cm -alpha60=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl5_alpha.nc'],'alpha_60cm',[j,i],[1,1]); -n60=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl5_n.nc'],'n_60cm',[j,i],[1,1]); -theta_s60=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl5_thetas.nc'],'thetas_60cm',[j,i],[1,1]); -theta_r60=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl5_thetar.nc'],'thetar_60cm',[j,i],[1,1]); -Ks60=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl5_Ks.nc'],'Ks_60cm',[j,i],[1,1]); -% 100cm -alpha100=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl6_alpha.nc'],'alpha_100cm',[j,i],[1,1]); -n100=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl6_n.nc'],'n_100cm',[j,i],[1,1]); -theta_s100=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl6_thetas.nc'],'thetas_100cm',[j,i],[1,1]); -theta_r100=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl6_thetar.nc'],'thetar_100cm',[j,i],[1,1]); -Ks100=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl6_Ks.nc'],'Ks_100cm',[j,i],[1,1]); -% 200cm -alpha200=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl7_alpha.nc'],'alpha_200cm',[j,i],[1,1]); -n200=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl7_n.nc'],'n_200cm',[j,i],[1,1]); -theta_s200=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl7_thetas.nc'],'thetas_200cm',[j,i],[1,1]); -theta_r200=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl7_thetar.nc'],'thetar_200cm',[j,i],[1,1]); -Ks200=ncread([SoilPropertyPath,'Schaap/PTF_SoilGrids_Schaap_sl7_Ks.nc'],'Ks_200cm',[j,i],[1,1]); - -%% load maximum fractional saturated area - -FMAX=ncread([SoilPropertyPath,'surfdata.nc'],'FMAX'); -i = round((latitude + 90)*2) + 1; - -if longitude < 0 - j = round((longitude + 360)*2) + 1; -else - j = round(longitude*2) + 1; -end -fmax=FMAX(j,i); - -% soil property -SaturatedK=[Ks0/(3600*24) Ks5/(3600*24) Ks30/(3600*24) Ks60/(3600*24) Ks100/(3600*24) Ks200/(3600*24)];%[2.67*1e-3 1.79*1e-3 1.14*1e-3 4.57*1e-4 2.72*1e-4]; %[2.3*1e-4 2.3*1e-4 0.94*1e-4 0.94*1e-4 0.68*1e-4] 0.18*1e-4Saturation hydraulic conductivity (cm.s^-1); -SaturatedMC=[theta_s0 theta_s5 theta_s30 theta_s60 theta_s100 theta_s200]; % 0.42 0.55 Saturated water content; -ResidualMC=[theta_r0 theta_r5 theta_r30 theta_r60 theta_r100 theta_r200]; %0.037 0.017 0.078 The residual water content of soil; -Coefficient_n=[n0 n5 n30 n60 n100 n200]; %1.2839 1.3519 1.2139 Coefficient in VG model; -Coefficient_Alpha=[alpha0 alpha5 alpha30 alpha60 alpha100 alpha200]; % 0.02958 0.007383 Coefficient in VG model; -porosity=[theta_s0 theta_s5 theta_s30 theta_s60 theta_s100 theta_s200]; % Soil porosity; -fieldMC=(1./(((341.09.*Coefficient_Alpha).^(Coefficient_n)+1).^(1-1./Coefficient_n))).*(SaturatedMC-ResidualMC)+ResidualMC; \ No newline at end of file diff --git a/utils/csv_to_nc/README.md b/utils/csv_to_nc/README.md index 2fb64964..487ecaf0 100644 --- a/utils/csv_to_nc/README.md +++ b/utils/csv_to_nc/README.md @@ -1,56 +1,11 @@ # Converting `.csv` files to NetCDF files -Currently, model outputs are several files in `csv` format. The model output -should be converted to one netcedf file according to Plumber protocol. To do so, -there is a file -[Variables_will_be_in_NetCDF_file.csv](./Variables_will_be_in_NetCDF_file.csv. -The file lists variables that should be in the netcdf file. Also, there is a -python script [csv_to_nc.py](./csv_to_nc.py) that contains main fucntions. Below -we explain how to use the python scripts. - -## 1. Create Conda environment - -> We need to do this step only once. - -We download and install Conda: - -```sh -wget https://github.com/conda-forge/miniforge/releases/latest/download/Mambaforge-pypy3-Linux-x86_64.sh -bash Mambaforge-pypy3-Linux-x86_64.sh -b -p ~/mamba -``` - -Then, update base environment: - -```sh -. ~/mamba/bin/activate -mamba update --name base mamba -``` - -Finally, we create new conda environment called 'stemmus' with all required dependencies: - -```sh -cd STEMMUS_SCOPE/utils/csv_to_nc -mamba env create -``` - -## 2. Activate Conda environment - -> We need to do this step before running our python scripts. - -The environment can be activated with - -```sh -. ~/mamba/bin/activate stemmus -``` - -## 3. Run python script - -Open the configuration file [config_file_crib.txt][../../config_file_crib.txt] -or [config_file_snellius.txt][../../config_file_snellius.txt] and edit paths. Then, - -```sh -cd STEMMUS_SCOPE -python utils/csv_to_nc/generate_netcdf_files.py --config_file config_file_crib.txt --variable_file utils/csv_to_nc/Variables_will_be_in_NetCDF_file.csv -``` - -This will generate `ECdata.csv` and a netcdf file related to model output. \ No newline at end of file +The model outputs are expected to be in one NetCDF file according to [ALMA +conventions](https://web.lmd.jussieu.fr/~polcher/ALMA/). The file +[required_netcdf_variables.csv](./required_netcdf_variables.csv) +lists required variable names and their attributes based on [`ALMA+CF` +convention table](https://docs.google.com/spreadsheets/d/1CA3aTvI9piXqRqO-3MGrsH1vW-Sd87D8iZXHGrqK42o/edit#gid=2085475627). + +The MATLAB source code writes model outputs in `csv` format in the output +directory. The NetCDF file is generated using the module `save.py` from +[`PyStemmusScope`](https://pystemmusscope.readthedocs.io/en/latest/autoapi/index.html). diff --git a/utils/csv_to_nc/__init__.py b/utils/csv_to_nc/__init__.py deleted file mode 100644 index e69de29b..00000000 diff --git a/utils/csv_to_nc/csv_to_nc.py b/utils/csv_to_nc/csv_to_nc.py deleted file mode 100644 index 189297eb..00000000 --- a/utils/csv_to_nc/csv_to_nc.py +++ /dev/null @@ -1,245 +0,0 @@ -import pandas as pd - -from netCDF4 import Dataset, num2date -from pathlib import Path - -def generate_ec_data(forcing_path, forcing_name, output_path, duration_size): - """ - Read forcing data and create EC data in csv format. - """ - - # read data - filename = Path(forcing_path, forcing_name) # care for windows/unix paths - nc_fid = Dataset(filename, mode='r') - - # Extract variables as numpy arrays - t = nc_fid.variables['time'][:].flatten()/3600/24 - - Ta = nc_fid.variables['Tair'][:].flatten()-273.15 - - Rin = nc_fid.variables['SWdown'][:].flatten() - - u = nc_fid.variables['Wind'][:].flatten() - - Rli = nc_fid.variables['LWdown'][:].flatten() - - p = nc_fid.variables['Psurf'][:].flatten()/100 - - RH = nc_fid.variables['RH'][:].flatten() - - LAI = nc_fid.variables['LAI'][:].flatten() - - nctime = nc_fid.variables['time'] - ncdate = num2date(nctime[:],units= nctime.units, calendar=nctime.calendar) - year = [date.year for date in ncdate] - - Pre = nc_fid.variables['Precip'][:].flatten()/10 - - CO2air = nc_fid.variables['CO2air'][:].flatten()*44/22.4 - - Qair = nc_fid.variables['Qair'][:].flatten() - - data = [t, Ta, Rin, u, Rli, p, RH, LAI, year, Pre, CO2air, Qair] - names = ['t', 'Ta', 'Rin', 'u', 'Rli', 'p', 'RH', 'LAI', 'year', 'Pre', 'CO2air', 'Qair'] - - # Write to csv - df = pd.DataFrame.from_dict(dict(zip(names, data))) - - # subset data, if needed - if duration_size != 'NA': - var_t_length = int(duration_size) - df = df.loc[:var_t_length-1] - - filename = Path(output_path, "ECdata.csv") - df.to_csv(filename, index=False) - - # Extra information - lat = float(nc_fid.variables['latitude'][:].flatten()[0]) - lon = float(nc_fid.variables['longitude'][:].flatten()[0]) - return filename, lat, lon, nctime - - -def split(s): - t=s.split('"') - u=[t[i] for i in range(0, len(t), 2)] - v=[t[i] for i in range(1, len(t), 2)] - #u=[u[i][(i%2):(len(u[i])-(i+1)%2)] for i in range(len(u))] - u=[ - u[i][(i%2):(len(u[i])-(i+1)%2)] for i in range(len(u)-1) - ] + [u[i][(i%2):] for i in range(len(u)-1, len(u))] - w = [] - for i in range(len(u)): - w.extend(u[i].split(',')) - if i < len(v): - w.append(v[i]) - return w - - -def readcsv(filename, nrHeaderLines): - f = open(filename) - header = f.readline() - header = header.strip().split(',') - if nrHeaderLines > 1: # it is either 1 or 2 - f.readline() - content = f.readlines() - data = {} - for line in content: - line = split(line.strip()) - for i in range(len(header)): - if header[i] != '': - if header[i] not in data: - data[header[i]] = [] - data[header[i]].append(line[i]) - return data - - -def read_depths(filename): - f = open(filename) - depths = f.readline() - depths = depths.strip().strip('#').strip(',').split() # the first line has ,,,,,, at the end - depths = [float(depth) for depth in depths] - return depths - - -def read2d_transposed_unit(filename, nrHeaderLines, unit, depths): - f = open(filename) - f.readline() # skip the headerline(s) - if nrHeaderLines > 1: # it is either 1 or 2 - f.readline() - content = f.readlines() - data = [] - for line in content: - line = line.strip().split(',') - line = [float(l) for l in line] # convert this to float as we may want to scale it - if unit == 'K': - # Celsius to Kelvin : K = 273.15 + C - line = [273.15 + c for c in line] - elif unit == 'kg/m2': - # Yijian Zeng: m3/m3 to kg/m2: SM = VolumetricWaterContent * Density * Thickness - # VolumetricWaterContent: provided (m3/m3) - # Density: constant (water_density = 1000 kg per m3) - # Thickness (m): compute from depth (cm) - line = [(1000.0 * vwc * depth / 100.0) for vwc,depth in zip(line, depths)] - data.append(line) - return data - - -def generateNetCdf(lat, lon, nctime, output_dir, csvfiles_path, variables_filename, duration_size): - # Renamed radiation.dat to radiation.csv - # Renamed LEtot to lEtot - # Split AResist into AResist_rac and AResist_ras - # Renamed the 2nd occurence of SWdown and LWdown to SWdown_ec and LWdown_ex - # Note that the values in this Excel sheet file determine the metadata that the variables will receive - - # specify additional metadata here: - - additional_metadata = { - 'model': 'STEMMUS_SCOPE', - 'institution': 'University of Twente; Northwest A&F University', - 'contact': 'Zhongbo Su, z.su@utwente.nl; Yijian Zeng, y.zeng@utwente.nl; Yunfei Wang, y.wang-3@utwente.nl', - 'license_type': 'CC BY 4.0', - 'license_url': 'https://creativecommons.org/licenses/by/4.0/', - 'latitude': lat, - 'longitude': lon - } - - # Our CSV reader can't guess the number of header-lines, so this is hardcoded here: - sim_theta = Path(csvfiles_path, 'Sim_Theta.csv') - sim_temp = Path(csvfiles_path, 'Sim_Temp.csv') - - headerlines = {'aerodyn.csv': 2, 'ECdata.csv': 1, 'fluxes.csv': 2, 'radiation.csv': 2, 'Sim_Temp.csv': 2, 'Sim_Theta.csv': 2, 'surftemp.csv': 2} - - print(f'Reading variable metadata from {variables_filename}') - variables = readcsv(variables_filename, 1) - depths = read_depths(sim_temp) - - # Create a new empty netCDF file, in NETCDF3_CLASSIC format, just like the example file AU-Tum_2002-2017_OzFlux_Met.nc - filename = f"{Path(csvfiles_path).stem}_STEMMUS_SCOPE.nc" - filename_out = Path(output_dir, filename) - print(f"Creating {filename_out} ") - nc = Dataset(filename_out, mode='w', format='NETCDF3_CLASSIC') - - # Create the dimensions, as required by netCDF - - nc.createDimension('x', size=1) - nc.createDimension('y', size=1) - nc.createDimension('z', size=len(depths)) - nc.createDimension('time', None) - nc.createDimension('nchar', size=200) # this is not used, however the example file has it - - # Create the variables, as required by netCDF - - var_x = nc.createVariable('x', 'float64', ('x')) - var_y = nc.createVariable('y', 'float64', ('y')) - var_z = nc.createVariable('z', 'float64', ('z')) - var_t = nc.createVariable('time', 'float64', ('time')) - var_latitude = nc.createVariable('latitude', 'float32', ('y','x')) - var_latitude.setncattr('long_name', 'Gridbox latitude') - var_latitude.setncattr('standard_name', 'latitude') - var_latitude.setncattr('units', 'degrees') - var_latitude[:] = lat - var_longitude = nc.createVariable('longitude', 'float32', ('y','x')) - var_longitude.setncattr('long_name', 'Gridbox longitude') - var_longitude.setncattr('standard_name', 'longitude') - var_longitude.setncattr('units', 'degrees') - var_longitude[:] = lon - - # Add the generic metadata (taken from additional_metadata above) - - for metadata in additional_metadata: - nc.setncattr(metadata, additional_metadata[metadata]) - - # Fill the x, y, time variables with values - - var_x[:] = lon # in the example AU-Tum_2002-2017_OzFlux_Met.nc this is 1 - var_y[:] = lat # in the example AU-Tum_2002-2017_OzFlux_Met.nc this is 1 - var_z.setncattr('units', 'depth in cm') - var_z[:] = depths - var_t.setncattr('units', f'{nctime.units}') - var_t.setncattr('calendar', f'{nctime.calendar}') - - # Add all parameters as a netCDF variable; also add the known metadata (units, long_name, Stemmus_name, definition) - - data = {} - - for i in range(len(variables['short_name_alma'])): - variable = variables['short_name_alma'][i] - file = variables['File name'][i] - stemmusname = variables['Variable name in STEMMUS-SCOPE'][i] - unit = variables['unit'][i] - long_name = variables['long_name'][i] - standard_name = variables['standard_name'][i] - definition = variables['definition'][i] - dimensions = variables['dimension'][i] - var = None - if dimensions == 'XYT': - var = nc.createVariable(variable, 'float32', ('time','y','x')) - elif dimensions == 'XYZT': - var = nc.createVariable(variable, 'float32', ('time','y','x','z')) - var.setncattr('units', unit) - var.setncattr('long_name', long_name) - var.setncattr('standard_name', standard_name) - if stemmusname != '': - var.setncattr('STEMMUS-SCOPE_name', stemmusname) - if definition != '': - var.setncattr('definition', definition) - if stemmusname != '': - if file not in data: - print(f'Reading data from file: {file}') - data[file] = readcsv(Path(csvfiles_path, file), headerlines[file]) - var[:] = data[file][stemmusname] - else: # Sim_Temp or Sim_Theta - print(f'Reading data from file: {file}') - matrix = read2d_transposed_unit(Path(csvfiles_path, file), headerlines[file], unit, depths) - var[:] = matrix - - # Finally fill var_t with the nr of seconds per timestep - # Note: we don't take the numbers from the file, because Year + DoY is not as accurate (it becomes 3599.99, 7199.99 etc) - if duration_size != 'NA': - var_t_length = int(duration_size) - else: - var_t_length = len(nctime[:]) - var_t[:] = [i*1800 for i in range(var_t_length)] - - nc.close() - return filename_out \ No newline at end of file diff --git a/utils/csv_to_nc/environment.yml b/utils/csv_to_nc/environment.yml deleted file mode 100644 index ec60be6a..00000000 --- a/utils/csv_to_nc/environment.yml +++ /dev/null @@ -1,8 +0,0 @@ ---- -name: stemmus -channels: - - conda-forge -dependencies: - - python=3.9 - - netcdf4>=1.5.8,<2 - - pandas diff --git a/utils/csv_to_nc/generate_netcdf_files.py b/utils/csv_to_nc/generate_netcdf_files.py deleted file mode 100644 index 6d016dc8..00000000 --- a/utils/csv_to_nc/generate_netcdf_files.py +++ /dev/null @@ -1,49 +0,0 @@ -"""Script to generate netcdf file from model output csv files. - -python generate_netcdf_files.py --config_file config_file_crib.txt --variable_file Variables_will_be_in_NetCDF_file.csv - -""" -import argparse -import csv_to_nc - -from pathlib import Path - - -def cli(): - """Generate netcdf file from model output csv files.""" - parser = argparse.ArgumentParser(description=__doc__) - parser.add_argument('--config_file', '-c', - default=Path(Path(__file__).parents[2],'config_file_crib.txt'), - help='Config file') - parser.add_argument('--variable_file', '-v', - default=Path(Path(__file__).parents[0],'Variables_will_be_in_NetCDF_file.csv'), - help='CSV file with variables that will be in NetCDF') - args = parser.parse_args() - - with open(args.config_file) as file: - config = dict(line.strip().split('=') for line in file) - - # get paths - output_path = config['OutputPath'] - forcing_path = config['ForcingPath'] - forcing_filename = config['ForcingFileName'] - - # get settings - duration_size = config['DurationSize'] - - # Find model output sub-directory linked to the forcing_filename - station_name = forcing_filename.split("_")[0] - model_output_dir = sorted(Path(output_path).glob(f"{station_name}*"))[-1] - - # Generate EC data from forcing data - ECdata_filename, lat, lon, nctime = csv_to_nc.generate_ec_data(forcing_path, forcing_filename, model_output_dir, duration_size) - print(f"{ECdata_filename}") - - # # Generate netcdf file from model output csv files - netcdf_filename = csv_to_nc.generateNetCdf(lat, lon, nctime, model_output_dir, model_output_dir, args.variable_file, duration_size) - print(f"Done writing {netcdf_filename}") - - -if __name__ == "__main__": - cli() - diff --git a/utils/csv_to_nc/Variables_will_be_in_NetCDF_file.csv b/utils/csv_to_nc/required_netcdf_variables.csv similarity index 97% rename from utils/csv_to_nc/Variables_will_be_in_NetCDF_file.csv rename to utils/csv_to_nc/required_netcdf_variables.csv index b26d3281..56950f25 100644 --- a/utils/csv_to_nc/Variables_will_be_in_NetCDF_file.csv +++ b/utils/csv_to_nc/required_netcdf_variables.csv @@ -1,27 +1,27 @@ -,pri_cmip,short_name_alma,short_name_cmip,standard_name,long_name,definition,unit,direction,dimension,grp_alma,grp_cmip,subgrid,Available in STEMMUS-SCOPE,File name,Variable name in STEMMUS-SCOPE, -1,1,SWnet,rss,surface_net_downward_shortwave_flux,Net shortwave radiation,"Incoming solar radiation less the simulated outgoing shortwave radiation, averaged over a grid cell",W/m2,Downward,XYT,,LEday,,Yes,radiation.csv,Netshort, -1,1,LWnet,rls,surface_net_downward_longwave_flux,Net longwave radiation,"Incident longwave radiation less the simulated outgoing longwave radiation, averaged over a grid cell",W/m2,Downward,XYT,,LEday,,Yes,radiation.csv,Netlong, -,2,SWdown,rsds,surface_downwelling_shortwave_flux_in_air,Downward short-wave radiation,,W/m2,Downward,XYT,,LEday,,Yes,radiation.csv,Rin, -,2,LWdown,rlds,surface_downwelling_longwave_flux_in_air,Downward long-wave radiation,,W/m2,Downward,XYT,,LEday,,Yes,radiation.csv,Rli, -,2,SWup,rsus,surface_upwelling_shortwave_flux_in_air,Upward short-wave radiation,,W/m2,Upward,XYT,,LEday,,Yes,radiation.csv,HemisOutShort, -2,2,LWup,rlus,surface_upwelling_longwave_flux_in_air,Upward long-wave radiation,This upward longwave flux is to be compared to an ISCCP derived product.,W/m2,Upward,XYT,,LEday,,Yes,radiation.csv,HemisOutLong, -1,1,Qle,hfls,surface_upward_latent_heat_flux,Latent heat flux,"Energy of evaporation, averaged over a grid cell",W/m2,Upward,XYT,,LEday,,Yes,fluxes.csv,lEtot, -1,1,Qh,hfss,surface_upward_sensible_heat_flux,Sensible heat flux,"Sensible energy, averaged over a grid cell",W/m2,Upward,XYT,,LEday,,Yes,fluxes.csv,Htot, -1,1,Qg,hfds,surface_downward_heat_flux,Ground heat flux,"Heat flux into the ground, averaged over a grid cell",W/m2,Downward,XYT,,LEday,,Yes,fluxes.csv,Gtot, -1,2,VegT,tcs,surface_canopy_skin_temperature,Vegetation Canopy Temperature,"Vegetation temperature, averaged over all vegetation types",K,-,XYT,,LEday,veg.,Yes,surftemp.csv,Tcave, -1,2,BaresoilT,tgs,surface_ground_skin_temperature,Temperature of bare soil,Surface bare soil temperature,K,-,XYT,,LEday,baresoil,Yes,surftemp.csv,Tsave, -2,1,SoilTemp,tsl,soil_temperature,Average layer soil temperature,Average soil temperature in each user-defined soil layer (3D variable),K,-,XYZT,,LEday,,Yes,Sim_Temp.csv,,"If soil layer thicknesses vary from one location to another, interpolate to a standard set of depths. Ideally, the interpolation should preserve the vertical integral." -1,1,SoilMoist,mrlsl,moisture_content_of_soil_layer,Average layer soil moisture,"Soil water content in each user-defined soil layer (3D variable). Includes the liquid, vapor and solid phases of water in the soil.",kg/m2,-,XYZT,,LWday,,Yes,Sim_Theta.csv,, -,2,AResist_rac,ares,aerodynamic_resistance,Aerodynamic resistance,,s/m,-,XYT,,LWday,,Yes,aerodyn.csv,rac, -,2,AResist_ras,ares,aerodynamic_resistance,Aerodynamic resistance,,s/m,-,XYT,,LWday,,Yes,aerodyn.csv,ras, -,1,RH,hur,relative_humidity,Relative humidity,,%,-,XYT,,LWday,,Yes,ECdata.csv,RH, -1,1,GPP,gpp,gross_primary_productivity_of_carbon,Gross Primary Production,Carbon Mass Flux out of Atmosphere due to Gross Primary Production on Land,kg/m2/s,Downward,XYT,,LCmon,,Yes,fluxes.csv,GPP, -1,1,SWdown_ec,rsds,surface_downwelling_shortwave_flux_in_air,Downward short-wave radiation,,W/m2,Downward,XYT,,L3hr,,,ECdata.csv,Rin, -1,1,LWdown_ec,rlds,surface_downwelling_longwave_flux_in_air,Downward long-wave radiation,,W/m2,Downward,XYT,,L3hr,,,ECdata.csv,Rli, -1,1,Qair,huss,specific_humidity,Near surface specific humidity,,kg/kg,-,XYT,,L3hr,,,ECdata.csv,Qair, -1,1,Tair,ta,air_temperature,Near surface air temperature,,K,-,XYT,,L3hr,,,ECdata.csv,Ta, -1,1,Psurf,ps,surface_air_pressure,Surface Pressure,,Pa,-,XYT,,L3hr,,,ECdata.csv,p, -2,1,Wind,ws,wind_speed,Near surface wind speed,,m/s,-,XYT,,L3hr,,,ECdata.csv,u, -,,Precip,pr,precipitation_flux,Precipitation rate,,kg/m2/s,Downward,XYT,,L3hr,,,ECdata.csv,Pre, -1,1,NEE,nep,surface_net_downward_mass_flux_of_carbon_dioxide_expressed_as_carbon_due_to_all_land_processes_excluding_anthropogenic_land_use_change,Net Ecosystem Exchange,Net Carbon Mass Flux out of Atmophere due to Net Ecosystem Productivity on Land.,kg/m2/s,Downward,XYT,,LCmon,,Yes,fluxes.csv,NEE, -1,1,Rnet,rss,surface_net_radiation_flux,Net radiation,"Net shortwave radiation less the simulated net longwave radiation, averaged over a grid cell",W/m2,Downward,XYT,,LEday,,Yes,radiation.csv,Rntot, +,pri_cmip,short_name_alma,short_name_cmip,standard_name,long_name,definition,unit,direction,dimension,grp_alma,grp_cmip,subgrid,Available in STEMMUS-SCOPE,file_name_STEMMUS-SCOPE,short_name_STEMMUS-SCOPE, +1,1,SWnet,rss,surface_net_downward_shortwave_flux,Net shortwave radiation,"Incoming solar radiation less the simulated outgoing shortwave radiation, averaged over a grid cell",W/m2,Downward,XYT,,LEday,,Yes,radiation.csv,Netshort, +1,1,LWnet,rls,surface_net_downward_longwave_flux,Net longwave radiation,"Incident longwave radiation less the simulated outgoing longwave radiation, averaged over a grid cell",W/m2,Downward,XYT,,LEday,,Yes,radiation.csv,Netlong, +,2,SWdown,rsds,surface_downwelling_shortwave_flux_in_air,Downward short-wave radiation,,W/m2,Downward,XYT,,LEday,,Yes,radiation.csv,Rin, +,2,LWdown,rlds,surface_downwelling_longwave_flux_in_air,Downward long-wave radiation,,W/m2,Downward,XYT,,LEday,,Yes,radiation.csv,Rli, +,2,SWup,rsus,surface_upwelling_shortwave_flux_in_air,Upward short-wave radiation,,W/m2,Upward,XYT,,LEday,,Yes,radiation.csv,HemisOutShort, +2,2,LWup,rlus,surface_upwelling_longwave_flux_in_air,Upward long-wave radiation,This upward longwave flux is to be compared to an ISCCP derived product.,W/m2,Upward,XYT,,LEday,,Yes,radiation.csv,HemisOutLong, +1,1,Qle,hfls,surface_upward_latent_heat_flux,Latent heat flux,"Energy of evaporation, averaged over a grid cell",W/m2,Upward,XYT,,LEday,,Yes,fluxes.csv,lEtot, +1,1,Qh,hfss,surface_upward_sensible_heat_flux,Sensible heat flux,"Sensible energy, averaged over a grid cell",W/m2,Upward,XYT,,LEday,,Yes,fluxes.csv,Htot, +1,1,Qg,hfds,surface_downward_heat_flux,Ground heat flux,"Heat flux into the ground, averaged over a grid cell",W/m2,Downward,XYT,,LEday,,Yes,fluxes.csv,Gtot, +1,2,VegT,tcs,surface_canopy_skin_temperature,Vegetation Canopy Temperature,"Vegetation temperature, averaged over all vegetation types",K,-,XYT,,LEday,veg.,Yes,surftemp.csv,Tcave, +1,2,BaresoilT,tgs,surface_ground_skin_temperature,Temperature of bare soil,Surface bare soil temperature,K,-,XYT,,LEday,baresoil,Yes,surftemp.csv,Tsave, +2,1,SoilTemp,tsl,soil_temperature,Average layer soil temperature,Average soil temperature in each user-defined soil layer (3D variable),K,-,XYZT,,LEday,,Yes,Sim_Temp.csv,,"If soil layer thicknesses vary from one location to another, interpolate to a standard set of depths. Ideally, the interpolation should preserve the vertical integral." +1,1,SoilMoist,mrlsl,moisture_content_of_soil_layer,Average layer soil moisture,"Soil water content in each user-defined soil layer (3D variable). Includes the liquid, vapor and solid phases of water in the soil.",kg/m2,-,XYZT,,LWday,,Yes,Sim_Theta.csv,, +,2,AResist_rac,ares,aerodynamic_resistance,Aerodynamic resistance,,s/m,-,XYT,,LWday,,Yes,aerodyn.csv,rac, +,2,AResist_ras,ares,aerodynamic_resistance,Aerodynamic resistance,,s/m,-,XYT,,LWday,,Yes,aerodyn.csv,ras, +,1,RH,hur,relative_humidity,Relative humidity,,%,-,XYT,,LWday,,Yes,ECdata.csv,RH, +1,1,GPP,gpp,gross_primary_productivity_of_carbon,Gross Primary Production,Carbon Mass Flux out of Atmosphere due to Gross Primary Production on Land,kg/m2/s,Downward,XYT,,LCmon,,Yes,fluxes.csv,GPP, +1,1,SWdown_ec,rsds,surface_downwelling_shortwave_flux_in_air,Downward short-wave radiation,,W/m2,Downward,XYT,,L3hr,,,ECdata.csv,Rin, +1,1,LWdown_ec,rlds,surface_downwelling_longwave_flux_in_air,Downward long-wave radiation,,W/m2,Downward,XYT,,L3hr,,,ECdata.csv,Rli, +1,1,Qair,huss,specific_humidity,Near surface specific humidity,,kg/kg,-,XYT,,L3hr,,,ECdata.csv,Qair, +1,1,Tair,ta,air_temperature,Near surface air temperature,,K,-,XYT,,L3hr,,,ECdata.csv,Ta, +1,1,Psurf,ps,surface_air_pressure,Surface Pressure,,Pa,-,XYT,,L3hr,,,ECdata.csv,p, +2,1,Wind,ws,wind_speed,Near surface wind speed,,m/s,-,XYT,,L3hr,,,ECdata.csv,u, +,,Precip,pr,precipitation_flux,Precipitation rate,,kg/m2/s,Downward,XYT,,L3hr,,,ECdata.csv,Pre, +1,1,NEE,nep,surface_net_downward_mass_flux_of_carbon_dioxide_expressed_as_carbon_due_to_all_land_processes_excluding_anthropogenic_land_use_change,Net Ecosystem Exchange,Net Carbon Mass Flux out of Atmophere due to Net Ecosystem Productivity on Land.,kg/m2/s,Downward,XYT,,LCmon,,Yes,fluxes.csv,NEE, +1,1,Rnet,rss,surface_net_radiation_flux,Net radiation,"Net shortwave radiation less the simulated net longwave radiation, averaged over a grid cell",W/m2,Downward,XYT,,LEday,,Yes,radiation.csv,Rntot,