diff --git a/.github/workflows/ci_run_scm_rts.yml b/.github/workflows/ci_run_scm_rts.yml index 632a44fbb..89db8a0a6 100644 --- a/.github/workflows/ci_run_scm_rts.yml +++ b/.github/workflows/ci_run_scm_rts.yml @@ -10,7 +10,7 @@ jobs: strategy: matrix: fortran-compiler: [gfortran-11] - build-type: [Release, Debug] + build-type: [Release, Debug, SinglePrecision] py-version: [3.9.12] # Environmental variables @@ -172,16 +172,30 @@ jobs: mkdir bin && cd bin cmake -DCCPP_SUITES=${suites},${suites_ps} -DCMAKE_BUILD_TYPE=Debug ../src + - name: Configure build with CMake (Single Precision) + if: matrix.build-type == 'SinglePrecision' + run: | + cd ${SCM_ROOT}/scm + mkdir bin && cd bin + cmake -DCCPP_SUITES=${suites},${suites_ps} -D32BIT=1 ../src + - name: Build SCM run: | cd ${SCM_ROOT}/scm/bin make -j4 - name: Run SCM RTs + if: matrix.build-type != 'SinglePrecision' run: | cd ${SCM_ROOT}/scm/bin ./run_scm.py --file /home/runner/work/ccpp-scm/ccpp-scm/test/rt_test_cases.py --runtime_mult 0.1 + - name: Run SCM Single Precision RTs + if: matrix.build-type == 'SinglePrecision' + run: | + cd ${SCM_ROOT}/scm/bin + ./run_scm.py --file /home/runner/work/ccpp-scm/ccpp-scm/test/rt_test_cases_sp.py --runtime_mult 0.1 + - name: Gather SCM RT output run: | cd ${SCM_ROOT}/test diff --git a/scm/doc/TechGuide/chap_quick.rst b/scm/doc/TechGuide/chap_quick.rst index 9f407dbba..f1524079e 100644 --- a/scm/doc/TechGuide/chap_quick.rst +++ b/scm/doc/TechGuide/chap_quick.rst @@ -5,7 +5,7 @@ Quick Start Guide This chapter provides instructions for obtaining and compiling the CCPP SCM. We provide instructions on building the code from scratch (:numref:`Section %s `), as well as -using Docker containers for machines that have Docker software installed (:numref:`Section %s `). +using Docker containers for machines that have Docker software installed (:numref:`Section %s `). .. _obtaining_code: @@ -15,7 +15,7 @@ Obtaining Code The source code for the SCM, CCPP, and their required components are provided through GitHub. The latest release branch contains the tested and supported version for general use, while the development branch (``main``) contains the latest -developer code, but may not be as stable or consistent with existing documentation. +developer code, but may not be as stable or consistent with existing documentation. Instructions for using either option are discussed here. Release Code @@ -29,7 +29,7 @@ Clone the source using The ``--recursive`` option is required to retrieve the ccpp-physics and ccpp-framework code, which are stored in separate repositories and linked to the SCM repository as submodules. -If not included initially, you can always retrieve the submodules +If not included initially, you can always retrieve the submodules by executing the following command from the SCM directory: .. code:: bash @@ -46,7 +46,7 @@ Development Code ^^^^^^^^^^^^^^^^ Developers seeking to contribute code to the SCM or CCPP will need to use the most up-to-date -version of the code, which can be found on the ``main`` branch of the repository: +version of the code, which can be found on the ``main`` branch of the repository: .. code:: bash @@ -121,7 +121,7 @@ System Requirements, Libraries, and Tools ----------------------------------------- The source code for the SCM and CCPP components is in the form of -programs written in FORTRAN 90 (with some required features from the +programs written in FORTRAN 90 (with some required features from the FORTRAN 2008 standard), and C. In addition, the model I/O relies on the NetCDF libraries, as well as the NCEP libraries ``bacio``, ``sp`` and ``w3emc``. @@ -156,7 +156,7 @@ contains the following set of libraries needed for building the SCM: Instructions for installing spack-stack can be found in the `spack-stack documentation `__. Spack-stack is already installed and maintained on many HPC platforms, including NSF NCAR's Derecho, NOAA's Hera and -Jet, and MSU's Orion. +Jet, and MSU's Orion. Compilers ^^^^^^^^^ @@ -230,7 +230,7 @@ For example, users on the NSF NCAR machine Derecho who wish to use Intel compile module load derecho_intel Additionally, for users who have installed spack-stack on their own MacOS or Linux machine can use the provided ``macos_clang`` -or ``linux_gnu`` modules. +or ``linux_gnu`` modules. .. note:: @@ -250,7 +250,7 @@ Python requirements The SCM build system invokes the ``ccpp_prebuild.py`` script, and so the Python environment must be set up prior to building. As mentioned earlier, a minimum Python version of 3.8 is required. Additionally, there are a few non-default modules required for the SCM to -function: ``f90nml`` (`documentation `__) and +function: ``f90nml`` (`documentation `__) and ``netcdf4`` (`documentation `__). Users can test if these are installed using this command in the shell: @@ -353,6 +353,13 @@ components. -DCMAKE_BUILD_TYPE=Debug + - Single Precision, lowers the default precision of variables from double to single precision. + A very few calculations are done in double precision where it is crucial to achieve results comparable to the default double precision. + + .. code:: bash + + -D32BIT=ON + - One can also save the output of this step to a log file: .. code:: bash @@ -425,7 +432,7 @@ Downloading input data The various SCM cases require staged input data in order to run. This includes input data for cases and lookup tables for runtime use. This is a large dataset (:math:`<`\ 1 GB) so it is not stored in the SCM repository, and must be downloaded -separately. To download this data place it in the correct directories, +separately. To download this data place it in the correct directories, execute the following scripts: .. code:: bash @@ -483,7 +490,7 @@ To see the full list of available options, use the ``--help`` flag: The run script’s full set of options are described below, where optional abbreviations are included in brackets. -If using the main branch, you should run the above command to ensure you have the most up-to-date list of options. +If using the main branch, you should run the above command to ensure you have the most up-to-date list of options. - ``--case [-c]`` @@ -750,7 +757,7 @@ Building the Docker image ^^^^^^^^^^^^^^^^^^^^^^^^^ The Dockerfile builds CCPP SCM v6.0.0 from source using the GNU -compiler. +compiler. The CCPP SCM has a number of system requirements and necessary libraries and tools. Below is a list, including versions, used to create the the @@ -794,7 +801,7 @@ and then executing the following steps: Inspect the Dockerfile if you would like to see details for how the image is built. The image will contain SCM prerequisite software from DTC, the SCM and CCPP code, and a pre-compiled executable for the SCM - with the 6 supported suites for the SCM. To view + with the 6 supported suites for the SCM. To view .. code:: bash @@ -874,7 +881,7 @@ Running the Docker image .. note:: Windows users may need to omit the curly braces around environment variables: use ``$OUT_DIR`` - instead of ``${OUT_DIR}``. + instead of ``${OUT_DIR}``. For running through all supported cases and suites, use @@ -914,7 +921,7 @@ Running the Docker image directory of the SCM with a pre-compiled executable. At this point, one could use the run scripts as described in previous sections (remembering to include the option on run scripts if output is to be - shared with the host machine). + shared with the host machine). .. warning:: diff --git a/scm/etc/scripts/precision_analysis.py b/scm/etc/scripts/precision_analysis.py new file mode 100644 index 000000000..2e4f8afd5 --- /dev/null +++ b/scm/etc/scripts/precision_analysis.py @@ -0,0 +1,338 @@ +import argparse +import itertools +import os +import subprocess +import sys +import xarray as xr +import numpy as np +import matplotlib.pyplot as plt + +# TODO: +# - [ ] add user argument to cmake "$@" in configure{32,64} +suites = ['SCM_RAP'] +cases = ['twpice'] + +global verbose +verbose = False + +def print_cmd(command): + print(f"$ {' '.join(command)}") + +def run_cmd(command): + print_cmd(command) + try: + # Run the cmake command + subprocess.run(command, check=True) + except subprocess.CalledProcessError as e: + # Handle errors if cmake command fails + print(f"Error running {e}") + sys.exit() + + +# Function to run configuration step +def configure32(): + print("Running single precision configuration step") + command = ["cmake", "src/", "-B", "build32", "-D32BIT=ON"] + run_cmd(command) + + +# Function to run configuration step +def configure64(): + print("Running double precision configuration step") + command = ["cmake", "src/", "-B", "build64"] + run_cmd(command) + + +# Function to run configuration step +def configure(): + print("Running configuration steps") + configure32() + configure64() + + +# Function to run build step +def build32(debug): + print("Running single precision build step") + if (debug): + command = ["make", "-C", "build32", "-j1", "VERBOSE=1"] + else: + command = ["make", "-C", "build32", "-j4"] + run_cmd(command) + +# Function to run build step +def build64(debug): + print("Running double precision build step") + if (debug): + command = ["make", "-C", "build64", "-j1", "VERBOSE=1"] + else: + command = ["make", "-C", "build64", "-j4"] + run_cmd(command) + + +# Function to run build step +def build(debug=False): + print("Running make steps") + build32(debug) + build64(debug) + + +# Post process single vs. double precision output +def post(): + print("Post processing single vs. double precision output") + postprocess() + + +def run32(): + print("Running single precision executable") + combinations = list(itertools.product(cases, suites)) + for (case, suite) in combinations: + command = ["build32/run_scm.py", "-c", case, + "-s", suite, + "--bin_dir", os.environ["PWD"]+"/build32", + "--run_dir", os.environ["PWD"]+"/output32"] + run_cmd(command) + +def run64(): + print("Running double precision executable") + combinations = list(itertools.product(cases, suites)) + for (case, suite) in combinations: + command = ["build64/run_scm.py", "-c", case, + "-s", suite, + "--bin_dir", os.environ["PWD"]+"/build64", + "--run_dir", os.environ["PWD"]+"/output64"] + run_cmd(command) + +# Function to run executable +def run(): + print("Running executables") + run32() + run64() + +def check_for_vertical_dimension(da): + vert_dim = None + if 'vert_dim_layer' in da.dims: + vert_dim = True + return vert_dim + + +def get_time_dimension(da): + # Initialize time variable to None + time = None + + # Iterate over dimensions of the Dataset + for dim in da.dims: + if 'time' in dim.lower(): + time = da[dim].name + break + + # Check if time variable was found + if time is None: + print(f"No 'time' dimension matched in {da.name}") + # else: + # print("Found 'time' dimension in the dataset.") + return time + +def check_dimensions(da1): + found_time = True + found_vert_layer = True + time_dim = get_time_dimension(da1) + vert_dim = check_for_vertical_dimension(da1) + if time_dim == None: + found_time = False + if verbose: + print(f"No time dimension found in variable {da1.name}") + return found_time + if vert_dim == None: + found_vert_layer = False + if verbose: + print(f"No vertical layer found in variable {da1.name}") + return found_vert_layer + return time_dim + +def calculate_mse(da1, da2, time_dim): + print("TIME DIM", time_dim) + """Calculate the mean square error (MSE) between two DataArrays.""" + mse = ((da1 - da2) ** 2)#.mean(dim=time_dim) + return mse + + +def save_mse_plots(case_dir, mse_dict): + """Save plots of mean square error (MSE) to a directory.""" + diff_dir = case_dir + '_diff' + os.makedirs(diff_dir, exist_ok=True) + + for var_name, mse in mse_dict.items(): + plt.figure() + mse.plot() + plt.title(f'Mean Square Error for {var_name}') + plt.xlabel('Column Level') + plt.ylabel('MSE') + plt.savefig(os.path.join(diff_dir, f'{var_name}_mse.png')) + plt.close() + +def save_re_plots(case_dir, re_dict): + """Save plots of relative error to a directory.""" + diff_dir = case_dir + '_diff' + os.makedirs(diff_dir, exist_ok=True) + + for var_name, re in re_dict.items(): + plt.figure() + re.plot() + plt.title(f'Relative Error for {var_name}') + plt.xlabel('Column Level') + plt.ylabel('Relative Error') + plt.savefig(os.path.join(diff_dir, f'{var_name}_re.png')) + plt.close() + + +def save_percentage_plots(case_dir, pd_dict): + """Save plots of percentage difference to a directory.""" + diff_dir = case_dir + '_diff' + os.makedirs(diff_dir, exist_ok=True) + + for var_name, pd in pd_dict.items(): + plt.figure() + pd.plot() + plt.title(f'Percentage Difference for {var_name}') + plt.xlabel('Column Level') + plt.ylabel('Percentage Diff (%)') + plt.savefig(os.path.join(diff_dir, f'{var_name}_pd.png')) + plt.close() + +# During post-processing step calculate relative the percentage difference +def calculate_relative_error_and_percentage(da32, da64, time_dim): + """Calculate the relative error between two DataArrays.""" + found_time = True + time_dim = get_time_dimension(da32) + if time_dim == None: + found_time = False + print(f"No time dimension found in variable {da32.name}") + return None, found_time + + relative_error = (np.abs(da64 - da32) / da64).mean(dim=time_dim) + percentage_dif = relative_error * 100 + return relative_error, percentage_dif + + +# Post processing function +def postprocess(): + print('Post-processing') + analysis_dir = 'output_diff' + analysis_file = 'variable_report.txt' + output_f = 'output.nc' + output32_dir = 'output32' + output64_dir = 'output64' + output32_dirs = sorted([os.path.basename(dirpath) for dirpath, _, _ in os.walk(output32_dir) if os.path.basename(dirpath).startswith('output_')]) + output64_dirs = sorted([os.path.basename(dirpath) for dirpath, _, _ in os.walk(output64_dir) if os.path.basename(dirpath).startswith('output_')]) + + vars_equal = [] + vars_diff = [] + + # Use set intersection to find directories that exist in both output32 and output64 + case_dirs = set(output32_dirs).intersection(output64_dirs) + for case_dir in case_dirs: + f32 = os.path.join(output32_dir, case_dir, output_f) + f64 = os.path.join(output64_dir, case_dir, output_f) + # Open f32 and f64 using xarray + ds32 = xr.open_dataset(f32) + ds64 = xr.open_dataset(f64) + + # Compare variables and compute mean square error (MSE) for variables that are different + mse_dict = {} + re_dict = {} + percentage_diff_dict = {} + for var_name in ds32.data_vars: + if var_name in ds64.data_vars: + if np.any(ds32[var_name] != ds64[var_name]): + time_dim = check_dimensions(ds32[var_name]) + if time_dim == False: + continue + mse = calculate_mse(ds32[var_name], ds64[var_name], time_dim) + relative_error, percentage_diff = \ + calculate_relative_error_and_percentage(ds32[var_name], + ds64[var_name], + time_dim) + mse_dict[var_name] = mse + re_dict[var_name] = relative_error + percentage_diff_dict[var_name] = percentage_diff + vars_diff.append(var_name) + diff_f = os.path.join(analysis_dir,case_dir) + print(f"Variable {var_name} differs in output files, writing diff to {diff_f}") + # Save plots of mean square error (MSE) to a directory + save_mse_plots(diff_f, mse_dict) + save_re_plots(diff_f, re_dict) + save_percentage_plots(diff_f, percentage_diff_dict) + else: + vars_equal.append(var_name) + print(f"Variable {var_name} is the same in both output files") + + # Write the vars_diff list to a text file + with open(os.path.join(analysis_dir, analysis_file), 'w') as f: + f.write("Different variables:\n") + for var in vars_diff: + f.write(f" - {var}\n") + f.write("Equal variables:\n") + for var in vars_equal: + f.write(f" - {var}\n") + + +def main(): + # Parse command-line arguments + parser = argparse.ArgumentParser(description="Script to configure, build, and run a program") + parser.add_argument("--configure", action="store_true", help="Configure CMake for single and double precision") + parser.add_argument("--configure32", action="store_true", help="Configure CMake for single precision") + parser.add_argument("--configure64", action="store_true", help="Configure CMake for double precision") + parser.add_argument("--build", action="store_true", help="Build single and double precision code") + parser.add_argument("--build32", action="store_true", help="Build single precision code") + parser.add_argument("--build64", action="store_true", help="Build double precision code") + parser.add_argument("--post", action="store_true", help="Postprocess the output") + parser.add_argument("--run", action="store_true", help="Run cases in single and double precision") + parser.add_argument("--run32", action="store_true", help="Run cases in single precision") + parser.add_argument("--run64", action="store_true", help="Run cases in double precision") + parser.add_argument("-d", "--debug", action="store_true", help="Debug mode") + parser.add_argument("-v", "--verbose", action="store_true", help="Verbose output") + args = parser.parse_args() + + # If --help is provided, print help message + if len(sys.argv) == 1 or "--help" in sys.argv: + parser.print_help() + sys.exit() + + verbose = args.verbose + + if args.debug: + debug = True + else: + debug = False + + # Execute requested step(s) + if args.configure: + configure() + elif args.configure32: + configure32() + elif args.configure64: + configure64() + elif args.build32: + build32(debug) + elif args.build64: + build64(debug) + elif args.build: + build(debug) + elif args.run: + run() + elif args.run32: + run32() + elif args.run64: + run64() + elif args.post: + post() + elif args.run: + run() + else: + configure() + build() + run() + + +if __name__ == "__main__": + main() diff --git a/scm/src/CMakeLists.txt b/scm/src/CMakeLists.txt index e8e8f7bbb..a42583c6c 100644 --- a/scm/src/CMakeLists.txt +++ b/scm/src/CMakeLists.txt @@ -50,7 +50,7 @@ else() RESULT_VARIABLE return_code ) endif() - + # Check return code from CCPP prebuild.py if(return_code EQUAL 0) message (STATUS "CCPP prebuild step completed successfully") @@ -71,6 +71,10 @@ find_package(NetCDF REQUIRED COMPONENTS C Fortran) find_package(bacio REQUIRED) find_package(sp REQUIRED) find_package(w3emc REQUIRED) +find_package(MPI REQUIRED) +if(NOT MPI_Fortran_HAVE_F08_MODULE) + message(FATAL_ERROR "MPI_F08 Required") +endif() SET(CCPP_FRAMEWORK_SRC ${CMAKE_SOURCE_DIR}/../../ccpp/framework) SET(CCPP_PHYSICS_SRC ${CMAKE_SOURCE_DIR}/../../ccpp/physics) @@ -90,6 +94,7 @@ endif() INCLUDE_DIRECTORIES(${CMAKE_BINARY_DIR}/ccpp/framework/src) INCLUDE_DIRECTORIES(${CMAKE_BINARY_DIR}/ccpp/physics) +include_directories(${MPI_Fortran_INCLUDE_PATH}) #------------------------------------------------------------------------------ # Add required preprocessor flags for build type @@ -131,17 +136,22 @@ set(SIMDMULTIARCH OFF CACHE BOOL "Enable multi-target SIMD instruction sets") #------------------------------------------------------------------------------ # Set compile options +if (32BIT) + add_definitions(-DSINGLE_PREC) + add_definitions(-DRTE_USE_SP) +endif() + if (${CMAKE_Fortran_COMPILER_ID} MATCHES "GNU") set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -ggdb -fbacktrace -cpp -fcray-pointer -ffree-line-length-none -fno-range-check") - + if(${CMAKE_Fortran_COMPILER_VERSION} VERSION_GREATER_EQUAL 10) set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fallow-argument-mismatch -fallow-invalid-boz") endif() - + if(NOT 32BIT) set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fdefault-real-8 -fdefault-double-8") endif() - + if (${CMAKE_BUILD_TYPE} MATCHES "Debug") set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -O0 -fno-unsafe-math-optimizations -frounding-math -fsignaling-nans -ffpe-trap=invalid,zero,overflow -fbounds-check") set(CMAKE_C_FLAGS_DEBUG "${CMAKE_C_FLAGS_DEBUG} -O0") @@ -149,7 +159,7 @@ if (${CMAKE_Fortran_COMPILER_ID} MATCHES "GNU") set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -O2") set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -O2") endif() - + set(CMAKE_C_FLAGS_RELEASE "-O2 -fPIC" CACHE STRING "" FORCE) set(CMAKE_Fortran_FLAGS_RELEASE "-O2 -fPIC" CACHE STRING "" FORCE) set(CMAKE_C_FLAGS_BITFORBIT "-O2 -fPIC" CACHE STRING "" FORCE) @@ -157,11 +167,11 @@ if (${CMAKE_Fortran_COMPILER_ID} MATCHES "GNU") elseif (${CMAKE_Fortran_COMPILER_ID} MATCHES "Intel") set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -g -traceback -fpp -fno-alias -auto -safe-cray-ptr -ftz -assume byterecl -nowarn -sox -align array64byte -qno-opt-dynamic-align") set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -qno-opt-dynamic-align -sox -fp-model source") - + if(NOT 32BIT) set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -real-size 64") endif() - + if (${CMAKE_BUILD_TYPE} MATCHES "Debug") set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -O0 -check -check noarg_temp_created -check nopointer -warn -warn noerrors -fp-stack-check -fstack-protector-all -fpe0 -debug -ftrapuv -init=snan,arrays") set(CMAKE_C_FLAGS_DEBUG "${CMAKE_C_FLAGS_DEBUG} -O0 -ftrapuv") @@ -184,7 +194,7 @@ elseif (${CMAKE_Fortran_COMPILER_ID} MATCHES "Intel") set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -march=core-avx-i") endif() endif() - + set(CMAKE_C_FLAGS_RELEASE "-O2 -fPIC" CACHE STRING "" FORCE) set(CMAKE_Fortran_FLAGS_RELEASE "-O2 -fPIC" CACHE STRING "" FORCE) set(CMAKE_C_FLAGS_BITFORBIT "-O2 -fPIC" CACHE STRING "" FORCE) @@ -202,7 +212,7 @@ elseif (${CMAKE_Fortran_COMPILER_ID} MATCHES "NVHPC") if(NOT 32BIT) set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -r8") endif() - + if (${CMAKE_BUILD_TYPE} MATCHES "Debug") set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -O0 -g") set(CMAKE_C_FLAGS_DEBUG "${CMAKE_C_FLAGS_DEBUG} -O0 -g") @@ -214,7 +224,7 @@ elseif (${CMAKE_Fortran_COMPILER_ID} MATCHES "NVHPC") set(MPI_C_COMPILER mpicc) set(MPI_CXX_COMPILER mpicxx) set(MPI_Fortran_COMPILER mpif90) - + set(CMAKE_C_FLAGS_RELEASE "-O2 -fPIC" CACHE STRING "" FORCE) set(CMAKE_Fortran_FLAGS_RELEASE "-O2 -fPIC" CACHE STRING "" FORCE) set(CMAKE_C_FLAGS_BITFORBIT "-O2 -fPIC" CACHE STRING "" FORCE) @@ -227,15 +237,17 @@ endif (${CMAKE_Fortran_COMPILER_ID} MATCHES "GNU") #------------------------------------------------------------------------------ # Set flag for 32bit dynamics build + if(32BIT) message(STATUS "Compile CCPP slow physics with 64-bit precision, fast physics with 32-bit precision") add_definitions(-DOVERLOAD_R4) if(CMAKE_Fortran_COMPILER_ID MATCHES "Intel") - set(CMAKE_Fortran_FLAGS_PHYSICS "-real-size 64 -no-prec-div -no-prec-sqrt") + # set(CMAKE_Fortran_FLAGS_PHYSICS "-real-size 64 -no-prec-div -no-prec-sqrt") elseif(CMAKE_Fortran_COMPILER_ID MATCHES "GNU") - set(CMAKE_Fortran_FLAGS_PHYSICS "-fdefault-real-8 -fdefault-double-8") + # set(CMAKE_Fortran_FLAGS_PHYSICS "-fdefault-real-8 -fdefault-double-8") endif() set(CMAKE_Fortran_FLAGS_DYNAMICS "") + else() message(STATUS "Compile CCPP physics with 64-bit precision") remove_definitions(-DOVERLOAD_R8) @@ -296,4 +308,3 @@ ADD_CUSTOM_COMMAND( COMMAND ${CMAKE_COMMAND} -E create_symlink ${CMAKE_SOURCE_DIR}/run_scm.py ${CMAKE_CURRENT_BINARY_DIR}/run_scm.py) - diff --git a/scm/src/scm_kinds.F90 b/scm/src/scm_kinds.F90 index 19fad9f25..1aad857e1 100644 --- a/scm/src/scm_kinds.F90 +++ b/scm/src/scm_kinds.F90 @@ -6,9 +6,17 @@ module scm_kinds !! \section arg_table_scm_kinds !! \htmlinclude scm_kinds.html !! - - integer, parameter :: sp = selected_real_kind(P= 6,R=37) - integer, parameter :: dp = selected_real_kind(P=13,R=300) + integer, parameter :: sp = selected_real_kind(P= 6,R=37) +#ifndef SINGLE_PREC + integer, parameter :: dp = selected_real_kind(P=13,R=300) integer, parameter :: qp = selected_real_kind(P=27,R=2400) +#else + integer, parameter :: dp = sp + integer, parameter :: qp = sp +#endif + ! these types exists to allow generic interface compilation in scm_utils.F90 + integer, parameter :: kind_scm_sp = selected_real_kind(P= 6,R=37) + integer, parameter :: kind_scm_dp = selected_real_kind(P=13,R=300) + integer, parameter :: kind_scm_qp = selected_real_kind(P=27,R=2400) end module scm_kinds diff --git a/scm/src/scm_utils.F90 b/scm/src/scm_utils.F90 index 48a33c4bf..ce2183afd 100644 --- a/scm/src/scm_utils.F90 +++ b/scm/src/scm_utils.F90 @@ -3,7 +3,7 @@ module scm_utils -use scm_kinds, only: sp, dp, qp +use scm_kinds, only: sp, dp, qp, kind_scm_dp use scm_physical_constants, only: con_rd, con_g implicit none @@ -87,9 +87,9 @@ subroutine find_vertical_index_pressure_sp(p_thresh, pres, k_out) real(kind=sp), intent(in) :: p_thresh real(kind=sp), intent(in) :: pres(:) integer, intent(out) :: k_out - + integer :: k - + k_out = -999 do k=1, size(pres) if (pres(k) <= p_thresh) then @@ -97,16 +97,16 @@ subroutine find_vertical_index_pressure_sp(p_thresh, pres, k_out) exit end if end do - + end subroutine find_vertical_index_pressure_sp subroutine find_vertical_index_pressure_dp(p_thresh, pres, k_out) - real(kind=dp), intent(in) :: p_thresh - real(kind=dp), intent(in) :: pres(:) + real(kind=kind_scm_dp), intent(in) :: p_thresh + real(kind=kind_scm_dp), intent(in) :: pres(:) integer, intent(out) :: k_out - + integer :: k - + k_out = -999 do k=1, size(pres) if (pres(k) <= p_thresh) then @@ -114,16 +114,16 @@ subroutine find_vertical_index_pressure_dp(p_thresh, pres, k_out) exit end if end do - + end subroutine find_vertical_index_pressure_dp subroutine find_vertical_index_height(z_thresh, height, k_out) real(kind=sp), intent(in) :: z_thresh real(kind=sp), intent(in) :: height(:) integer, intent(out) :: k_out - + integer k - + k_out = -999 do k=1, size(height) if (height(k) >= z_thresh) then @@ -131,7 +131,7 @@ subroutine find_vertical_index_height(z_thresh, height, k_out) exit end if end do - + end subroutine find_vertical_index_height subroutine w_to_omega(n_col, n_lev, w, p, T, omega) @@ -153,7 +153,7 @@ integer function lcm(a,b) integer, intent(in) :: a,b lcm = a*b / gcd(a,b) end function lcm - + integer function gcd(a,b) integer, intent(in) :: a,b integer :: c,d,t @@ -172,15 +172,15 @@ end function gcd end module scm_utils module NetCDF_read - use scm_kinds, only : sp, dp, qp + use scm_kinds, only : sp, dp, qp, kind_scm_dp use netcdf - + implicit none - + real(kind=sp) :: missing_value = -9999.0 integer :: missing_value_int = -9999 character (3) :: missing_value_char = "mis" - + interface NetCDF_read_var module procedure NetCDF_read_var_0d_int module procedure NetCDF_read_var_1d_int @@ -197,32 +197,32 @@ module NetCDF_read module procedure NetCDF_read_var_3d_dp module procedure NetCDF_read_var_4d_dp end interface - + interface NetCDF_conditionally_read_var module procedure NetCDF_conditionally_read_var_1d_sp module procedure NetCDF_conditionally_read_var_2d_sp module procedure NetCDF_conditionally_read_var_3d_sp module procedure NetCDF_conditionally_read_var_4d_sp end interface NetCDF_conditionally_read_var - + interface NetCDF_read_att module procedure NetCDF_read_att_char module procedure NetCDF_read_att_int module procedure NetCDF_read_att_char_or_int module procedure NetCDF_read_att_sp end interface - + contains - + subroutine NetCDF_read_var_0d_int(ncid, var_name, req, var_data) - + integer, intent(in) :: ncid character (*), intent(in) :: var_name logical, intent(in) :: req integer, intent(out) :: var_data - + integer :: varID, ierr - + if (req) then call check(NF90_INQ_VARID(ncid,var_name,varID),var_name) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) @@ -234,18 +234,18 @@ subroutine NetCDF_read_var_0d_int(ncid, var_name, req, var_data) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) end if end if - + end subroutine NetCDF_read_var_0d_int - + subroutine NetCDF_read_var_1d_int(ncid, var_name, req, var_data) - + integer, intent(in) :: ncid character (*), intent(in) :: var_name logical, intent(in) :: req integer, dimension(:), intent(out) :: var_data - + integer :: varID, ierr - + if (req) then call check(NF90_INQ_VARID(ncid,var_name,varID),var_name) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) @@ -257,18 +257,18 @@ subroutine NetCDF_read_var_1d_int(ncid, var_name, req, var_data) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) end if end if - + end subroutine NetCDF_read_var_1d_int - + subroutine NetCDF_read_var_2d_int(ncid, var_name, req, var_data) - + integer, intent(in) :: ncid character (*), intent(in) :: var_name logical, intent(in) :: req integer, dimension(:,:), intent(out) :: var_data - + integer :: varID, ierr - + if (req) then call check(NF90_INQ_VARID(ncid,var_name,varID),var_name) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) @@ -280,18 +280,18 @@ subroutine NetCDF_read_var_2d_int(ncid, var_name, req, var_data) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) end if end if - + end subroutine NetCDF_read_var_2d_int - + subroutine NetCDF_read_var_3d_int(ncid, var_name, req, var_data) - + integer, intent(in) :: ncid character (*), intent(in) :: var_name logical, intent(in) :: req integer, dimension(:,:,:), intent(out) :: var_data - + integer :: varID, ierr - + if (req) then call check(NF90_INQ_VARID(ncid,var_name,varID),var_name) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) @@ -303,18 +303,18 @@ subroutine NetCDF_read_var_3d_int(ncid, var_name, req, var_data) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) end if end if - + end subroutine NetCDF_read_var_3d_int - + subroutine NetCDF_read_var_0d_sp(ncid, var_name, req, var_data) - + integer, intent(in) :: ncid character (*), intent(in) :: var_name logical, intent(in) :: req real(kind=sp), intent(out) :: var_data - + integer :: varID, ierr - + if (req) then call check(NF90_INQ_VARID(ncid,var_name,varID),var_name) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) @@ -326,18 +326,18 @@ subroutine NetCDF_read_var_0d_sp(ncid, var_name, req, var_data) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) end if end if - + end subroutine NetCDF_read_var_0d_sp - + subroutine NetCDF_read_var_1d_sp(ncid, var_name, req, var_data) - + integer, intent(in) :: ncid character (*), intent(in) :: var_name logical, intent(in) :: req real(kind=sp), dimension(:), intent(out) :: var_data - + integer :: varID, ierr - + if (req) then call check(NF90_INQ_VARID(ncid,var_name,varID),var_name) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) @@ -349,18 +349,18 @@ subroutine NetCDF_read_var_1d_sp(ncid, var_name, req, var_data) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) end if end if - + end subroutine NetCDF_read_var_1d_sp - + subroutine NetCDF_read_var_2d_sp(ncid, var_name, req, var_data) - + integer, intent(in) :: ncid character (*), intent(in) :: var_name logical, intent(in) :: req real(kind=sp), dimension(:,:), intent(out) :: var_data - + integer :: varID, ierr - + if (req) then call check(NF90_INQ_VARID(ncid,var_name,varID),var_name) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) @@ -372,18 +372,18 @@ subroutine NetCDF_read_var_2d_sp(ncid, var_name, req, var_data) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) end if end if - + end subroutine NetCDF_read_var_2d_sp - + subroutine NetCDF_read_var_3d_sp(ncid, var_name, req, var_data) - + integer, intent(in) :: ncid character (*), intent(in) :: var_name logical, intent(in) :: req real(kind=sp), dimension(:,:,:), intent(out) :: var_data - + integer :: varID, ierr - + if (req) then call check(NF90_INQ_VARID(ncid,var_name,varID),var_name) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) @@ -395,18 +395,18 @@ subroutine NetCDF_read_var_3d_sp(ncid, var_name, req, var_data) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) end if end if - + end subroutine NetCDF_read_var_3d_sp - + subroutine NetCDF_read_var_4d_sp(ncid, var_name, req, var_data) - + integer, intent(in) :: ncid character (*), intent(in) :: var_name logical, intent(in) :: req real(kind=sp), dimension(:,:,:,:), intent(out) :: var_data - + integer :: varID, ierr - + if (req) then call check(NF90_INQ_VARID(ncid,var_name,varID),var_name) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) @@ -418,18 +418,18 @@ subroutine NetCDF_read_var_4d_sp(ncid, var_name, req, var_data) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) end if end if - + end subroutine NetCDF_read_var_4d_sp - + subroutine NetCDF_read_var_0d_dp(ncid, var_name, req, var_data) - + integer, intent(in) :: ncid character (*), intent(in) :: var_name logical, intent(in) :: req - real(kind=dp), intent(out) :: var_data - + real(kind=kind_scm_dp), intent(out) :: var_data + integer :: varID, ierr - + if (req) then call check(NF90_INQ_VARID(ncid,var_name,varID),var_name) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) @@ -441,18 +441,18 @@ subroutine NetCDF_read_var_0d_dp(ncid, var_name, req, var_data) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) end if end if - + end subroutine NetCDF_read_var_0d_dp - + subroutine NetCDF_read_var_1d_dp(ncid, var_name, req, var_data) - + integer, intent(in) :: ncid character (*), intent(in) :: var_name logical, intent(in) :: req - real(kind=dp), dimension(:), intent(out) :: var_data - + real(kind=kind_scm_dp), dimension(:), intent(out) :: var_data + integer :: varID, ierr - + if (req) then call check(NF90_INQ_VARID(ncid,var_name,varID),var_name) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) @@ -464,18 +464,18 @@ subroutine NetCDF_read_var_1d_dp(ncid, var_name, req, var_data) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) end if end if - + end subroutine NetCDF_read_var_1d_dp - + subroutine NetCDF_read_var_2d_dp(ncid, var_name, req, var_data) - + integer, intent(in) :: ncid character (*), intent(in) :: var_name logical, intent(in) :: req - real(kind=dp), dimension(:,:), intent(out) :: var_data - + real(kind=kind_scm_dp), dimension(:,:), intent(out) :: var_data + integer :: varID, ierr - + if (req) then call check(NF90_INQ_VARID(ncid,var_name,varID),var_name) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) @@ -487,18 +487,18 @@ subroutine NetCDF_read_var_2d_dp(ncid, var_name, req, var_data) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) end if end if - + end subroutine NetCDF_read_var_2d_dp - + subroutine NetCDF_read_var_3d_dp(ncid, var_name, req, var_data) - + integer, intent(in) :: ncid character (*), intent(in) :: var_name logical, intent(in) :: req - real(kind=dp), dimension(:,:,:), intent(out) :: var_data - + real(kind=kind_scm_dp), dimension(:,:,:), intent(out) :: var_data + integer :: varID, ierr - + if (req) then call check(NF90_INQ_VARID(ncid,var_name,varID),var_name) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) @@ -510,18 +510,18 @@ subroutine NetCDF_read_var_3d_dp(ncid, var_name, req, var_data) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) end if end if - + end subroutine NetCDF_read_var_3d_dp - + subroutine NetCDF_read_var_4d_dp(ncid, var_name, req, var_data) - + integer, intent(in) :: ncid character (*), intent(in) :: var_name logical, intent(in) :: req - real(kind=dp), dimension(:,:,:,:), intent(out) :: var_data - + real(kind=kind_scm_dp), dimension(:,:,:,:), intent(out) :: var_data + integer :: varID, ierr - + if (req) then call check(NF90_INQ_VARID(ncid,var_name,varID),var_name) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) @@ -533,9 +533,9 @@ subroutine NetCDF_read_var_4d_dp(ncid, var_name, req, var_data) call check(NF90_GET_VAR(ncid,varID,var_data),var_name) end if end if - + end subroutine NetCDF_read_var_4d_dp - + !Generic subroutine to check for netCDF I/O errors subroutine check(status, var_name) integer, intent ( in) :: status @@ -546,16 +546,16 @@ subroutine check(status, var_name) stop "stopped" end if end subroutine check - + subroutine NetCDF_read_att_int(ncid, var_id, att_name, req, att_data) - + integer, intent(in) :: ncid, var_id character (*), intent(in) :: att_name logical, intent(in) :: req integer, intent(out) :: att_data - + integer :: varID, ierr - + if (req) then ierr = NF90_INQUIRE_ATTRIBUTE(ncid, var_id, att_name) if (ierr /= NF90_NOERR) then @@ -572,18 +572,18 @@ subroutine NetCDF_read_att_int(ncid, var_id, att_name, req, att_data) call check(NF90_GET_ATT(ncid, var_id, att_name, att_data),att_name) end if end if - + end subroutine NetCDF_read_att_int - + subroutine NetCDF_read_att_sp(ncid, var_id, att_name, req, att_data) - + integer, intent(in) :: ncid, var_id character (*), intent(in) :: att_name logical, intent(in) :: req real(kind=sp), intent(out) :: att_data - + integer :: varID, ierr - + if (req) then ierr = NF90_INQUIRE_ATTRIBUTE(ncid, var_id, att_name) if (ierr /= NF90_NOERR) then @@ -600,18 +600,18 @@ subroutine NetCDF_read_att_sp(ncid, var_id, att_name, req, att_data) call check(NF90_GET_ATT(ncid, var_id, att_name, att_data),att_name) end if end if - + end subroutine NetCDF_read_att_sp - + subroutine NetCDF_read_att_char(ncid, var_id, att_name, req, att_data) - + integer, intent(in) :: ncid, var_id character (*), intent(in) :: att_name logical, intent(in) :: req character (*), intent(out) :: att_data - + integer :: varID, ierr - + if (req) then ierr = NF90_INQUIRE_ATTRIBUTE(ncid, var_id, att_name) if (ierr /= NF90_NOERR) then @@ -628,19 +628,19 @@ subroutine NetCDF_read_att_char(ncid, var_id, att_name, req, att_data) call check(NF90_GET_ATT(ncid, var_id, att_name, att_data),att_name) end if end if - + end subroutine NetCDF_read_att_char - + subroutine NetCDF_read_att_char_or_int(ncid, var_id, att_name, req, att_data, char_att_data) - + integer, intent(in) :: ncid, var_id character (*), intent(in) :: att_name logical, intent(in) :: req character (*), intent(out) :: char_att_data integer, intent(out) :: att_data - + integer :: varID, ierr, type - + if (req) then ierr = NF90_INQUIRE_ATTRIBUTE(ncid, var_id, att_name, xtype = type) if (ierr /= NF90_NOERR) then @@ -676,7 +676,7 @@ subroutine NetCDF_read_att_char_or_int(ncid, var_id, att_name, req, att_data, ch end if end if end if - + end subroutine NetCDF_read_att_char_or_int subroutine NetCDF_conditionally_read_var_1d_sp(var_ctl, var_att, var_name, filename, ncid, var_data) @@ -715,16 +715,16 @@ subroutine NetCDF_conditionally_read_var_2d_sp(var_ctl, var_att, var_name, filen else var_data = missing_value end if - end subroutine NetCDF_conditionally_read_var_2d_sp + end subroutine NetCDF_conditionally_read_var_2d_sp subroutine NetCDF_conditionally_read_var_3d_sp(var_ctl, var_att, var_name, filename, ncid, var_data) integer, intent(in) :: var_ctl, ncid character (*), intent(in) :: var_att, var_name, filename real(kind=sp), dimension(:,:,:), intent(out) :: var_data real(kind=sp) :: missing_value_eps - + missing_value_eps = missing_value + 0.01 - + if (var_ctl > 0) then call NetCDF_read_var(ncid, var_name, .False., var_data) if (maxval(var_data) < missing_value_eps) then @@ -735,15 +735,15 @@ subroutine NetCDF_conditionally_read_var_3d_sp(var_ctl, var_att, var_name, filen var_data = missing_value end if end subroutine NetCDF_conditionally_read_var_3d_sp - + subroutine NetCDF_conditionally_read_var_4d_sp(var_ctl, var_att, var_name, filename, ncid, var_data) integer, intent(in) :: var_ctl, ncid character (*), intent(in) :: var_att, var_name, filename real(kind=sp), dimension(:,:,:,:), intent(out) :: var_data real(kind=sp) :: missing_value_eps - + missing_value_eps = missing_value + 0.01 - + if (var_ctl > 0) then call NetCDF_read_var(ncid, var_name, .False., var_data) if (maxval(var_data) < missing_value_eps) then @@ -754,17 +754,17 @@ subroutine NetCDF_conditionally_read_var_4d_sp(var_ctl, var_att, var_name, filen var_data = missing_value end if end subroutine NetCDF_conditionally_read_var_4d_sp - -end module NetCDF_read + +end module NetCDF_read module NetCDF_def use NetCDF_read, only : check use netcdf - + implicit none - + contains - + subroutine NetCDF_def_var(ncid, var_name, var_type, desc, unit, varid, dims) use NetCDF_read, only: missing_value, missing_value_int integer, intent(in) :: ncid @@ -774,7 +774,7 @@ subroutine NetCDF_def_var(ncid, var_name, var_type, desc, unit, varid, dims) character (*), intent(in) :: desc character (*), intent(in) :: unit integer, intent(out) :: varid - + if (present(dims)) then call CHECK(NF90_DEF_VAR(NCID=ncid,NAME=var_name,XTYPE=var_type,DIMIDS=dims,VARID=varid),var_name) else @@ -782,7 +782,7 @@ subroutine NetCDF_def_var(ncid, var_name, var_type, desc, unit, varid, dims) end if call CHECK(NF90_PUT_ATT(NCID=ncid,VARID=varid,NAME="description",VALUES=desc),var_name) call CHECK(NF90_PUT_ATT(NCID=ncid,VARID=varid,NAME="units",VALUES=unit),var_name) - + if (var_type == NF90_FLOAT) then call CHECK(NF90_PUT_ATT(NCID=ncid,VARID=varid,NAME="_FillValue",VALUES=missing_value),var_name) elseif (var_type == NF90_INT) then @@ -791,7 +791,7 @@ subroutine NetCDF_def_var(ncid, var_name, var_type, desc, unit, varid, dims) write(0,'(a,i0,a)') "The variable '" // var_name // "' is defined as a type other than NF90_FLOAT or NF90_INT. Stopping..." STOP end if - + end subroutine NetCDF_def_var end module NetCDF_def @@ -799,25 +799,25 @@ module NetCDF_put use NetCDF_read, only : check use netcdf use scm_kinds, only : sp, dp, qp - + implicit none - + interface NetCDF_put_var module procedure NetCDF_put_var_int_0d module procedure NetCDF_put_var_1d module procedure NetCDF_put_var_2d end interface - + contains - + subroutine NetCDF_put_var_int_0d(ncid, var_name, var, known_varid) integer, intent(in) :: ncid character (*), intent(in) :: var_name integer, intent(in) :: var integer, intent(in), optional :: known_varid - + integer :: var_id - + !write(*,*) 'Putting variable: ',var_name if (present(known_varid)) then CALL CHECK(NF90_PUT_VAR(NCID=ncid,VARID=known_varid,VALUES=var),var_name) @@ -825,17 +825,17 @@ subroutine NetCDF_put_var_int_0d(ncid, var_name, var, known_varid) CALL CHECK(NF90_INQ_VARID(NCID=ncid,NAME=var_name,VARID=var_id),var_name) CALL CHECK(NF90_PUT_VAR(NCID=ncid,VARID=var_id,VALUES=var),var_name) end if - + end subroutine NetCDF_put_var_int_0d - + subroutine NetCDF_put_var_1d(ncid, var_name, var, itt, mult_const) integer, intent(in) :: ncid, itt character (*), intent(in) :: var_name real(kind=dp), intent(in), dimension(:) :: var real(kind=dp), intent(in), optional :: mult_const - + integer :: var_id - + !write(*,*) 'Putting variable: ',var_name CALL CHECK(NF90_INQ_VARID(NCID=ncid,NAME=var_name,VARID=var_id),var_name) if (present(mult_const)) then @@ -843,17 +843,17 @@ subroutine NetCDF_put_var_1d(ncid, var_name, var, itt, mult_const) else CALL CHECK(NF90_PUT_VAR(NCID=ncid,VARID=var_id,VALUES=var,START=(/1,itt /)),var_name) end if - + end subroutine NetCDF_put_var_1d - + subroutine NetCDF_put_var_2d(ncid, var_name, var, itt, mult_const) integer, intent(in) :: ncid, itt character (*), intent(in) :: var_name real(kind=dp), intent(in), dimension(:,:) :: var real(kind=dp), intent(in), optional :: mult_const - + integer :: var_id - + !write(*,*) 'Putting variable: ',var_name CALL CHECK(NF90_INQ_VARID(NCID=ncid,NAME=var_name,VARID=var_id),var_name) if (present(mult_const)) then @@ -861,73 +861,73 @@ subroutine NetCDF_put_var_2d(ncid, var_name, var, itt, mult_const) else CALL CHECK(NF90_PUT_VAR(NCID=ncid,VARID=var_id,VALUES=var,START=(/1,1,itt /)),var_name) end if - + end subroutine NetCDF_put_var_2d end module NetCDF_put module data_qc use scm_kinds, only : sp, dp, qp use NetCDF_read, only: missing_value, missing_value_int - + implicit none - + interface check_missing module procedure check_missing_int_0d module procedure check_missing_int_1d module procedure check_missing_real_dp_0d module procedure check_missing_real_dp_1d end interface - + interface conditionally_set_var module procedure conditionally_set_int_var_0d module procedure conditionally_set_int_var_1d module procedure conditionally_set_real_dp_var_0d module procedure conditionally_set_real_dp_var_1d end interface - + contains - + subroutine check_missing_int_0d(var, missing) integer, intent(in) :: var logical, intent(out) :: missing - + missing = .false. if (var == missing_value_int) missing = .true. end subroutine check_missing_int_0d - + subroutine check_missing_int_1d(var, missing) integer, dimension(:), intent(in) :: var logical, intent(out) :: missing - + missing = .false. if ( ANY(var == missing_value_int)) missing = .true. end subroutine check_missing_int_1d - + subroutine check_missing_real_dp_0d(var, missing) real(kind=dp), intent(in) :: var logical, intent(out) :: missing - + missing = .false. if (var == missing_value) missing = .true. end subroutine check_missing_real_dp_0d - + subroutine check_missing_real_dp_1d(var, missing) real(kind=dp), dimension(:), intent(in) :: var logical, intent(out) :: missing - + missing = .false. if ( ANY(var == missing_value)) missing = .true. end subroutine check_missing_real_dp_1d - + subroutine conditionally_set_int_var_0d(input, set_var, input_name, req, missing) integer, intent(in) :: input integer, intent(inout) :: set_var character (*), intent(in) :: input_name logical, intent(in) :: req logical, intent(out) :: missing - + call check_missing(input, missing) - + if (.not. missing) then set_var = input else @@ -936,18 +936,18 @@ subroutine conditionally_set_int_var_0d(input, set_var, input_name, req, missing STOP end if end if - + end subroutine conditionally_set_int_var_0d - + subroutine conditionally_set_int_var_1d(input, set_var, input_name, req, missing) integer, dimension(:), intent(in) :: input integer, dimension(:), intent(inout) :: set_var character (*), intent(in) :: input_name logical, intent(in) :: req logical, intent(out) :: missing - + call check_missing(input, missing) - + if (.not. missing) then set_var = input else @@ -956,18 +956,18 @@ subroutine conditionally_set_int_var_1d(input, set_var, input_name, req, missing STOP end if end if - + end subroutine conditionally_set_int_var_1d - + subroutine conditionally_set_real_dp_var_0d(input, set_var, input_name, req, missing) real(kind=dp), intent(in) :: input real(kind=dp), intent(inout) :: set_var character (*), intent(in) :: input_name logical, intent(in) :: req logical, intent(out) :: missing - + call check_missing(input, missing) - + if (.not. missing) then set_var = input else @@ -976,18 +976,18 @@ subroutine conditionally_set_real_dp_var_0d(input, set_var, input_name, req, mis STOP end if end if - + end subroutine conditionally_set_real_dp_var_0d - + subroutine conditionally_set_real_dp_var_1d(input, set_var, input_name, req, missing) real(kind=dp), dimension(:), intent(in) :: input real(kind=dp), dimension(:), intent(inout) :: set_var character (*), intent(in) :: input_name logical, intent(in) :: req logical, intent(out) :: missing - + call check_missing(input, missing) - + if (.not. missing) then set_var = input else @@ -996,6 +996,6 @@ subroutine conditionally_set_real_dp_var_1d(input, set_var, input_name, req, mis STOP end if end if - + end subroutine conditionally_set_real_dp_var_1d end module data_qc diff --git a/scm/src/scm_vgrid.F90 b/scm/src/scm_vgrid.F90 index f827f89f9..12ac24fc0 100644 --- a/scm/src/scm_vgrid.F90 +++ b/scm/src/scm_vgrid.F90 @@ -3,7 +3,7 @@ module scm_vgrid -use scm_kinds, only: sp, dp, qp +use scm_kinds, only: sp, dp, qp, kind_scm_dp, kind_scm_sp use scm_physical_constants, only : con_cp, con_rocp, con_fvirt, con_g, con_rd implicit none @@ -48,6 +48,13 @@ subroutine get_FV3_vgrid(scm_input, scm_state) character(len=80) :: line character(len=16) :: file_format + integer :: nx, ny + real(kind_scm_dp), allocatable :: pres_l_row(:), pres_i(:,:) + real(kind_scm_dp), parameter :: zero_dp = 0.0 + ! added for forcing initialized pressure to be single precision for + ! single and double precision runs + real(kind_scm_sp), parameter :: zero_sp = 0.0 + real(kind_scm_sp), allocatable :: pres_l_row_sp(:) #include "fv_eta.h" @@ -539,17 +546,25 @@ subroutine get_FV3_vgrid(scm_input, scm_state) p_ref = scm_input%input_pres_surf(1) pres_sfc_inv = 1.0/p_ref + + nx = size(scm_state%pres_i, 1) + ny = size(scm_state%pres_i, 2) + allocate(pres_i(nx,ny), source=zero_dp) do k=1, km+1 - scm_state%pres_i(:,k) = scm_state%a_k(k) + scm_state%b_k(k)*p_ref + pres_i(:,k) = scm_state%a_k(k) + scm_state%b_k(k)*p_ref scm_state%si(:,k) = scm_state%a_k(k)*pres_sfc_inv + scm_state%b_k(k) scm_state%exner_i(:,k) = (scm_state%pres_i(:,k)/1.0E5)**con_rocp end do + scm_state%pres_i = pres_i !> - Calculate layer center pressures, sigma, and exner function. + allocate(pres_l_row(nx), source=zero_dp) + allocate(pres_l_row_sp(nx), source=zero_sp) do k=1, km - scm_state%pres_l(:,k) = ((1.0/(con_rocp+1.0))*& - (scm_state%pres_i(:,k)**(con_rocp+1.0) - scm_state%pres_i(:,k+1)**(con_rocp+1.0))/ & - (scm_state%pres_i(:,k) - scm_state%pres_i(:,k+1)))**(1.0/con_rocp) + pres_l_row_sp = ((1.0/(con_rocp+1.0))*& + (pres_i(:,k)**(con_rocp+1.0) - pres_i(:,k+1)**(con_rocp+1.0))/ & + (pres_i(:,k) - pres_i(:,k+1)))**(1.0/con_rocp) + scm_state%pres_l(:,k) = pres_l_row_sp scm_state%sl(:,k) = 0.5*(scm_state%si(:,k) + scm_state%si(:,k+1)) scm_state%exner_l(:,k) = (scm_state%pres_l(:,k)/1.0E5)**con_rocp diff --git a/test/rt_test_cases_sp.py b/test/rt_test_cases_sp.py new file mode 100644 index 000000000..a9a888642 --- /dev/null +++ b/test/rt_test_cases_sp.py @@ -0,0 +1,10 @@ +run_list = [\ + #---------------------------------------------------------------------------------------------------------------------------------------------- + # CCPP-SCM single precision supported suites + #---------------------------------------------------------------------------------------------------------------------------------------------- + {"case": "arm_sgp_summer_1997_A", "suite": "SCM_RAP"}, \ + {"case": "twpice", "suite": "SCM_RAP"}, \ + {"case": "bomex", "suite": "SCM_RAP"}, \ + {"case": "astex", "suite": "SCM_RAP"}, \ + {"case": "LASSO_2016051812", "suite": "SCM_RAP"}, \ + ]