From 5a91ce60aebc1d8c026cc7c1fa3d75e958c09f9a Mon Sep 17 00:00:00 2001 From: Suvarchal Kumar Cheedela Date: Tue, 23 Jun 2026 03:04:25 +0200 Subject: [PATCH 1/8] add single-precision (WP=4) build support for standalone fesom.x Make working precision configurable via -DUSE_SINGLE_PRECISION=ON. Default build stays double precision and unchanged. - WP is now preprocessor-selected (4/8) in o_PARAM; central MPI_WP parameter in MOD_PARTIT mirrors it (MPI_REAL vs MPI_DOUBLE_PRECISION). - CMake gains USE_SINGLE_PRECISION option with per-compiler default-real flag flip (Intel -r4, GNU drops -fdefault-real-8, Cray -s real32, NVHPC -Mr4); CVMix forced OFF in SP (double-only, out of scope). - Replace MPI_DOUBLE_PRECISION/MPI_DOUBLE with MPI_WP at WP-data sites; convert WP comm/gather buffers from real(real64) to real(kind=WP). - NetCDF stays double on disk; WP reads routed through generic nf_*_vara_x so existing meshes/restarts/forcing keep working (convert-on-I/O). - mod_transit and reference-density locals converted to WP; io_meandata vector rotation round-trips through WP temps; iceberg subsystem made WP-consistent. - Startup banner prints WP byte size and SINGLE/DOUBLE mode. Co-authored-by: Cursor --- CMakeLists.txt | 7 ++ src/CMakeLists.txt | 47 +++++++++++--- src/MOD_PARTIT.F90 | 6 ++ src/cavity_param.F90 | 6 +- src/fesom_module.F90 | 5 ++ src/gen_halo_exchange.F90 | 106 +++++++++++++++---------------- src/gen_ic3d.F90 | 72 ++++++++++++--------- src/gen_modules_partitioning.F90 | 24 +++---- src/gen_modules_read_NetCDF.F90 | 42 ++++++------ src/gen_support.F90 | 12 ++-- src/gen_surface_forcing.F90 | 18 +++--- src/icb_allocate.F90 | 2 +- src/icb_elem.F90 | 22 +++---- src/icb_modules.F90 | 2 +- src/icb_step.F90 | 16 ++--- src/icb_thermo.F90 | 78 +++++++++++------------ src/io_fesom_file.F90 | 20 +++--- src/io_meandata.F90 | 9 ++- src/io_restart_file_group.F90 | 16 ++--- src/mod_transit.F90 | 53 ++++++++-------- src/oce_ale_pressure_bv.F90 | 10 +-- src/oce_dyn.F90 | 4 +- src/oce_mesh.F90 | 28 ++++---- src/oce_modules.F90 | 7 +- src/solver.F90 | 8 +-- src/write_step_info.F90 | 82 ++++++++++++------------ 26 files changed, 385 insertions(+), 317 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 0d4ca2784..b7b5b6414 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -27,6 +27,13 @@ set(VERBOSE OFF CACHE BOOL "toggle debug output") set(FESOM_PROFILING OFF CACHE BOOL "Enable enhanced profiling system") set(CVMIX ON CACHE BOOL "Download/Compile/Install/Link CVMix libray") +# Single precision (WP=4) build: CVMix is a double-precision library and is out of scope +# for the SP port, so disable it to keep the SP build clean (native KPP/PP still work). +if(USE_SINGLE_PRECISION) + message(STATUS "USE_SINGLE_PRECISION=ON: forcing CVMIX=OFF (out of scope for single precision)") + set(CVMIX OFF CACHE BOOL "Download/Compile/Install/Link CVMix libray" FORCE) +endif() + # Testing options option(BUILD_TESTING "Build tests" OFF) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c4cdbe4f2..d36a3392d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -79,6 +79,9 @@ message(STATUS "RECOM_COUPLED: ${RECOM_COUPLED}") option(USE_ICEPACK "Use ICEPACK" OFF) message(STATUS "USE_ICEPACK: ${USE_ICEPACK}") +option(USE_SINGLE_PRECISION "Build FESOM in single precision (WP=4)" OFF) +message(STATUS "USE_SINGLE_PRECISION: ${USE_SINGLE_PRECISION}") + # option to trigger building a library version of FESOM # we do not always build the library along with the executable to avoid having two targets here in the CMakeLists.txt @@ -292,6 +295,10 @@ if(USE_ICEPACK) target_compile_definitions(${PROJECT_NAME} PRIVATE __icepack) endif() +if(USE_SINGLE_PRECISION) + target_compile_definitions(${PROJECT_NAME} PRIVATE USE_SINGLE_PRECISION) +endif() + if(CVMIX) target_compile_definitions(${PROJECT_NAME} PRIVATE __cvmix) endif() @@ -323,8 +330,12 @@ endif() # CMAKE_Fortran_COMPILER_ID will also work if a wrapper is being used (e.g. mpif90 wraps ifort -> compiler id is Intel) if(${CMAKE_Fortran_COMPILER_ID} STREQUAL Intel OR ${CMAKE_Fortran_COMPILER_ID} STREQUAL IntelLLVM ) - # Base compiler flags - target_compile_options(${PROJECT_NAME} PRIVATE -O3 -r8 -i4 -fp-model precise -no-prec-div -fimf-use-svml -init=zero -no-wrap-margin -fpe0 -fpp) + # Base compiler flags (default real kind follows USE_SINGLE_PRECISION) + if(USE_SINGLE_PRECISION) + target_compile_options(${PROJECT_NAME} PRIVATE -O3 -r4 -i4 -fp-model precise -no-prec-div -fimf-use-svml -init=zero -no-wrap-margin -fpe0 -fpp) + else() + target_compile_options(${PROJECT_NAME} PRIVATE -O3 -r8 -i4 -fp-model precise -no-prec-div -fimf-use-svml -init=zero -no-wrap-margin -fpe0 -fpp) + endif() # compiler flags not supported by IntelLLVM if(${CMAKE_Fortran_COMPILER_ID} STREQUAL Intel ) @@ -381,11 +392,17 @@ if(${CMAKE_Fortran_COMPILER_ID} STREQUAL Intel OR ${CMAKE_Fortran_COMPILER_ID} elseif(${CMAKE_Fortran_COMPILER_ID} STREQUAL GNU ) # target_compile_options(${PROJECT_NAME} PRIVATE -O3 -finit-local-zero -finline-functions -fimplicit-none -fdefault-real-8 -ffree-line-length-none) + # default real kind follows USE_SINGLE_PRECISION: promote to 8 for DP, leave at 4 for SP + if(USE_SINGLE_PRECISION) + set(_FESOM_GNU_RKIND "") + else() + set(_FESOM_GNU_RKIND -fdefault-real-8 -fdefault-double-8) + endif() if(${FESOM_PLATFORM_STRATEGY} STREQUAL ubuntu ) message(STATUS "Allowing type mismatches on Ubuntu for CI Testing" ) # NOTE(PG): Would be nicer to grab the CI=True from the env variable - target_compile_options(${PROJECT_NAME} PRIVATE -O2 -g -fbacktrace -ffloat-store -finit-local-zero -finline-functions -fimplicit-none -fdefault-real-8 -fdefault-double-8 -ffree-line-length-none -cpp) + target_compile_options(${PROJECT_NAME} PRIVATE -O2 -g -fbacktrace -ffloat-store -finit-local-zero -finline-functions -fimplicit-none ${_FESOM_GNU_RKIND} -ffree-line-length-none -cpp) else() - target_compile_options(${PROJECT_NAME} PRIVATE -O3 -ffloat-store -finit-local-zero -finline-functions -fimplicit-none -fdefault-real-8 -fdefault-double-8 -ffree-line-length-none -cpp) + target_compile_options(${PROJECT_NAME} PRIVATE -O3 -ffloat-store -finit-local-zero -finline-functions -fimplicit-none ${_FESOM_GNU_RKIND} -ffree-line-length-none -cpp) endif() if(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") if(CMAKE_Fortran_COMPILER_VERSION VERSION_GREATER_EQUAL 10) @@ -450,7 +467,13 @@ elseif(${CMAKE_Fortran_COMPILER_ID} STREQUAL GNU ) elseif(${CMAKE_Fortran_COMPILER_ID} STREQUAL Cray ) #target_compile_options(${PROJECT_NAME} PRIVATE -c -emf -hbyteswapio -hflex_mp=conservative -hfp1 -hadd_paren -Ounroll0 -hipa0 -r am -s real64 -N 1023 -g -G2 -O3) - target_compile_options(${PROJECT_NAME} PRIVATE -c -emf -hbyteswapio -hflex_mp=conservative -hfp1 -hadd_paren -Ounroll0 -hipa0 -r am -s real64 -N 1023 -g -G2 -O2 -hnoacc -M878) #-hnoacc is a workaround for cray automatically activate -hacc, -M878 is to suppress ftn-878 warning + # default real kind follows USE_SINGLE_PRECISION + if(USE_SINGLE_PRECISION) + set(_FESOM_CRAY_RKIND real32) + else() + set(_FESOM_CRAY_RKIND real64) + endif() + target_compile_options(${PROJECT_NAME} PRIVATE -c -emf -hbyteswapio -hflex_mp=conservative -hfp1 -hadd_paren -Ounroll0 -hipa0 -r am -s ${_FESOM_CRAY_RKIND} -N 1023 -g -G2 -O2 -hnoacc -M878) #-hnoacc is a workaround for cray automatically activate -hacc, -M878 is to suppress ftn-878 warning if(${ENABLE_OPENMP}) target_compile_options(${PROJECT_NAME} PRIVATE -homp) else() @@ -464,10 +487,16 @@ elseif(${CMAKE_Fortran_COMPILER_ID} STREQUAL Cray ) endif() elseif(${CMAKE_Fortran_COMPILER_ID} STREQUAL NVHPC ) target_compile_definitions(${PROJECT_NAME} PRIVATE ENABLE_NVHPC_WORKAROUNDS) + # default real kind follows USE_SINGLE_PRECISION + if(USE_SINGLE_PRECISION) + set(_FESOM_NVHPC_RKIND -Mr4) + else() + set(_FESOM_NVHPC_RKIND -Mr8) + endif() if(${ENABLE_OPENACC}) target_compile_options(${PROJECT_NAME} PRIVATE - $<$:-Mallocatable=95 -Mr8 -pgf90libs -Mnofma -Minfo=all -acc=verystrict -gpu=math_uniform,cuda12.2,cc80> - $<$:-Mallocatable=95 -Mr8 -pgf90libs -Minfo=all -acc=verystrict -gpu=cc80>) + $<$:-Mallocatable=95 ${_FESOM_NVHPC_RKIND} -pgf90libs -Mnofma -Minfo=all -acc=verystrict -gpu=math_uniform,cuda12.2,cc80> + $<$:-Mallocatable=95 ${_FESOM_NVHPC_RKIND} -pgf90libs -Minfo=all -acc=verystrict -gpu=cc80>) set(CMAKE_EXE_LINKER_FLAGS_DEBUG "-acc=verystrict -Mnofma -gpu=math_uniform,cuda12.2,cc80") set(CMAKE_EXE_LINKER_FLAGS_RELEASE "-acc=verystrict -gpu=cc80") if(${DISABLE_OPENACC_ATOMICS}) @@ -475,8 +504,8 @@ elseif(${CMAKE_Fortran_COMPILER_ID} STREQUAL NVHPC ) endif() else() target_compile_options(${PROJECT_NAME} PRIVATE - $<$:-Mallocatable=95 -Mr8 -pgf90libs -Mnofma> - $<$:-Mallocatable=95 -Mr8 -pgf90libs>) + $<$:-Mallocatable=95 ${_FESOM_NVHPC_RKIND} -pgf90libs -Mnofma> + $<$:-Mallocatable=95 ${_FESOM_NVHPC_RKIND} -pgf90libs>) endif() endif() diff --git a/src/MOD_PARTIT.F90 b/src/MOD_PARTIT.F90 index 2e8330a4b..3c4d81f28 100644 --- a/src/MOD_PARTIT.F90 +++ b/src/MOD_PARTIT.F90 @@ -11,6 +11,12 @@ module MOD_PARTIT #endif IMPLICIT NONE SAVE +! MPI datatype matching the working precision WP (see o_PARAM) +#if defined(USE_SINGLE_PRECISION) +integer, parameter :: MPI_WP = MPI_REAL ! single precision +#else +integer, parameter :: MPI_WP = MPI_DOUBLE_PRECISION ! double precision (default) +#endif integer, parameter :: MAX_LAENDERECK=16 integer, parameter :: MAX_NEIGHBOR_PARTITIONS=32 diff --git a/src/cavity_param.F90 b/src/cavity_param.F90 index 1c2a42132..b022e3a1a 100644 --- a/src/cavity_param.F90 +++ b/src/cavity_param.F90 @@ -131,9 +131,9 @@ subroutine compute_nrst_pnt2cavline(partit, mesh) ! mpi reduce from local to global call MPI_AllREDUCE(lcl_cavl_idx, cavl_idx, nod2d, MPI_INTEGER , MPI_SUM, MPI_COMM_FESOM, MPIerr) - call MPI_AllREDUCE(lcl_cavl_lon, cavl_lon, nod2d, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) - call MPI_AllREDUCE(lcl_cavl_lat, cavl_lat, nod2d, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) - call MPI_AllREDUCE(lcl_cavl_dep, cavl_dep, nod2d, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(lcl_cavl_lon, cavl_lon, nod2d, MPI_WP, MPI_SUM, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(lcl_cavl_lat, cavl_lat, nod2d, MPI_WP, MPI_SUM, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(lcl_cavl_dep, cavl_dep, nod2d, MPI_WP, MPI_SUM, MPI_COMM_FESOM, MPIerr) !___________________________________________________________________________ ! deallocate local arrays diff --git a/src/fesom_module.F90 b/src/fesom_module.F90 index 9ea20678a..52cb049bb 100755 --- a/src/fesom_module.F90 +++ b/src/fesom_module.F90 @@ -212,6 +212,11 @@ subroutine fesom_init(fesom_total_nsteps) print *,"FESOM2 git SHA: "//fesom_git_sha() call MPI_Get_library_version(f%mpi_version_txt, f%mpi_version_len, f%MPIERR) print *,"MPI library version: "//trim(f%mpi_version_txt) +#if defined(USE_SINGLE_PRECISION) + print '(a,i0,a)'," FESOM working precision: WP=",WP," bytes (SINGLE PRECISION MODE)" +#else + print '(a,i0,a)'," FESOM working precision: WP=",WP," bytes (DOUBLE PRECISION MODE)" +#endif print *, achar(27)//'[32m' //'____________________________________________________________'//achar(27)//'[0m' print *, achar(27)//'[7;32m'//' --> FESOM BUILDS UP MODEL CONFIGURATION '//achar(27)//'[0m' end if diff --git a/src/gen_halo_exchange.F90 b/src/gen_halo_exchange.F90 index 92b9aba02..1afe6108c 100755 --- a/src/gen_halo_exchange.F90 +++ b/src/gen_halo_exchange.F90 @@ -146,7 +146,7 @@ subroutine exchange_nod2D(nod_array2D, partit, luse_g2g) USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod_array2D(:) +real(kind=WP), intent(inout) :: nod_array2D(:) logical, intent(in),optional :: luse_g2g #include "associate_part_def.h" #include "associate_part_ass.h" @@ -166,7 +166,7 @@ subroutine exchange_nod2D_begin(nod_array2D, partit, luse_g2g) USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod_array2D(:) +real(kind=WP), intent(inout) :: nod_array2D(:) integer :: n, sn, rn logical, intent(in),optional :: luse_g2g logical :: lg2g @@ -227,8 +227,8 @@ subroutine exchange_nod2D_2fields(nod1_array2D, nod2_array2D, partit, luse_g2g) USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod1_array2D(:) -real(real64), intent(inout) :: nod2_array2D(:) +real(kind=WP), intent(inout) :: nod1_array2D(:) +real(kind=WP), intent(inout) :: nod2_array2D(:) logical, intent(in),optional :: luse_g2g #include "associate_part_def.h" #include "associate_part_ass.h" @@ -248,8 +248,8 @@ subroutine exchange_nod2D_2fields_begin(nod1_array2D, nod2_array2D, partit, luse USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod1_array2D(:) -real(real64), intent(inout) :: nod2_array2D(:) +real(kind=WP), intent(inout) :: nod1_array2D(:) +real(kind=WP), intent(inout) :: nod2_array2D(:) integer :: n, sn, rn logical, intent(in),optional :: luse_g2g logical :: lg2g @@ -322,9 +322,9 @@ subroutine exchange_nod2D_3fields(nod1_array2D, nod2_array2D, nod3_array2D, part USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod1_array2D(:) -real(real64), intent(inout) :: nod2_array2D(:) -real(real64), intent(inout) :: nod3_array2D(:) +real(kind=WP), intent(inout) :: nod1_array2D(:) +real(kind=WP), intent(inout) :: nod2_array2D(:) +real(kind=WP), intent(inout) :: nod3_array2D(:) logical, intent(in),optional :: luse_g2g #include "associate_part_def.h" #include "associate_part_ass.h" @@ -345,9 +345,9 @@ subroutine exchange_nod2D_3fields_begin(nod1_array2D, nod2_array2D, nod3_array2D USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod1_array2D(:) -real(real64), intent(inout) :: nod2_array2D(:) -real(real64), intent(inout) :: nod3_array2D(:) +real(kind=WP), intent(inout) :: nod1_array2D(:) +real(kind=WP), intent(inout) :: nod2_array2D(:) +real(kind=WP), intent(inout) :: nod3_array2D(:) integer :: n, sn, rn logical, intent(in),optional :: luse_g2g logical :: lg2g @@ -435,7 +435,7 @@ subroutine exchange_nod3D(nod_array3D, partit, luse_g2g) USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod_array3D(:,:) +real(kind=WP), intent(inout) :: nod_array3D(:,:) logical, intent(in),optional :: luse_g2g if (partit%npes > 1) then @@ -454,7 +454,7 @@ subroutine exchange_nod3D_begin(nod_array3D, partit, luse_g2g) USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod_array3D(:,:) +real(kind=WP), intent(inout) :: nod_array3D(:,:) integer :: n, sn, rn integer :: nz, nl1 logical, intent(in),optional :: luse_g2g @@ -525,8 +525,8 @@ subroutine exchange_nod3D_2fields(nod1_array3D,nod2_array3D, partit, luse_g2g) USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod1_array3D(:,:) -real(real64), intent(inout) :: nod2_array3D(:,:) +real(kind=WP), intent(inout) :: nod1_array3D(:,:) +real(kind=WP), intent(inout) :: nod2_array3D(:,:) logical, intent(in),optional :: luse_g2g #include "associate_part_def.h" #include "associate_part_ass.h" @@ -546,8 +546,8 @@ subroutine exchange_nod3D_2fields_begin(nod1_array3D,nod2_array3D, partit, luse_ USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod1_array3D(:,:) -real(real64), intent(inout) :: nod2_array3D(:,:) +real(kind=WP), intent(inout) :: nod1_array3D(:,:) +real(kind=WP), intent(inout) :: nod2_array3D(:,:) integer :: n, sn, rn integer :: nz, nl1, nl2 logical, intent(in),optional :: luse_g2g @@ -638,7 +638,7 @@ subroutine exchange_nod3D_n(nod_array3D, partit, luse_g2g) USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod_array3D(:,:,:) +real(kind=WP), intent(inout) :: nod_array3D(:,:,:) logical, intent(in),optional :: luse_g2g if (partit%npes>1) then call exchange_nod3D_n_begin(nod_array3D, partit, luse_g2g) @@ -656,7 +656,7 @@ subroutine exchange_nod3D_n_begin(nod_array3D, partit, luse_g2g) USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: nod_array3D(:,:,:) +real(kind=WP), intent(inout) :: nod_array3D(:,:,:) integer :: n, sn, rn integer :: nz, nl1, n_val logical, intent(in),optional :: luse_g2g @@ -772,7 +772,7 @@ subroutine exchange_elem3D(elem_array3D, partit, luse_g2g) USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: elem_array3D(:,:) +real(kind=WP), intent(inout) :: elem_array3D(:,:) logical, intent(in),optional :: luse_g2g #include "associate_part_def.h" #include "associate_part_ass.h" @@ -790,7 +790,7 @@ subroutine exchange_elem3D_begin(elem_array3D, partit, luse_g2g) USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: elem_array3D(:,:) +real(kind=WP), intent(inout) :: elem_array3D(:,:) integer :: n, sn, rn, nl1 logical :: lg2g logical, intent(in),optional :: luse_g2g @@ -991,7 +991,7 @@ subroutine exchange_elem3D_n(elem_array3D, partit, luse_g2g) USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: elem_array3D(:,:,:) +real(kind=WP), intent(inout) :: elem_array3D(:,:,:) logical, intent(in),optional :: luse_g2g #include "associate_part_def.h" #include "associate_part_ass.h" @@ -1010,7 +1010,7 @@ subroutine exchange_elem3D_n_begin(elem_array3D, partit, luse_g2g) USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: elem_array3D(:,:,:) +real(kind=WP), intent(inout) :: elem_array3D(:,:,:) integer :: n, sn, rn, n_val, nl1 logical :: lg2g logical, intent(in),optional :: luse_g2g @@ -1129,7 +1129,7 @@ subroutine exchange_elem2D(elem_array2D, partit, luse_g2g) USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: elem_array2D(:) +real(kind=WP), intent(inout) :: elem_array2D(:) logical, intent(in),optional :: luse_g2g #include "associate_part_def.h" #include "associate_part_ass.h" @@ -1149,7 +1149,7 @@ subroutine exchange_elem2D_begin(elem_array2D, partit, luse_g2g) USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(inout), target :: partit -real(real64), intent(inout) :: elem_array2D(:) +real(kind=WP), intent(inout) :: elem_array2D(:) integer :: n, sn, rn logical :: lg2g logical, intent(in),optional :: luse_g2g @@ -1338,11 +1338,11 @@ subroutine broadcast_nod3D(arr3D, arr3Dglobal, partit) type(t_partit), intent(inout), target :: partit INTEGER :: nz, counter,nl1 integer :: i, n, nTS, sender, status(MPI_STATUS_SIZE) -real(real64) :: arr3D(:,:) -real(real64) :: arr3Dglobal(:,:) +real(kind=WP) :: arr3D(:,:) +real(kind=WP) :: arr3Dglobal(:,:) integer :: node_size INTEGER, ALLOCATABLE, DIMENSION(:) :: irecvbuf -real(real64), ALLOCATABLE, DIMENSION(:) :: sendbuf, recvbuf +real(kind=WP), ALLOCATABLE, DIMENSION(:) :: sendbuf, recvbuf #include "associate_part_def.h" #include "associate_part_ass.h" @@ -1368,7 +1368,7 @@ subroutine broadcast_nod3D(arr3D, arr3Dglobal, partit) ENDDO ENDDO - CALL MPI_SEND(sendbuf(1), nTS*nl1, MPI_DOUBLE_PRECISION, & + CALL MPI_SEND(sendbuf(1), nTS*nl1, MPI_WP, & sender, 2, MPI_COMM_FESOM, MPIerr ) DEALLOCATE(irecvbuf, sendbuf) @@ -1379,7 +1379,7 @@ subroutine broadcast_nod3D(arr3D, arr3Dglobal, partit) MPI_COMM_FESOM, MPIerr ) ALLOCATE(recvbuf(node_size*nl1)) - CALL MPI_RECV( recvbuf(1), node_size*nl1, MPI_DOUBLE_PRECISION, 0, & + CALL MPI_RECV( recvbuf(1), node_size*nl1, MPI_WP, 0, & 2, MPI_COMM_FESOM, status, MPIerr ) counter=0 DO n = 1, node_size @@ -1403,11 +1403,11 @@ subroutine broadcast_nod2D(arr2D, arr2Dglobal, partit) USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(in), target :: partit -real(real64) :: arr2D(:) -real(real64) :: arr2Dglobal(:) +real(kind=WP) :: arr2D(:) +real(kind=WP) :: arr2Dglobal(:) integer :: i, n, nTS, sender, status(MPI_STATUS_SIZE) INTEGER, ALLOCATABLE, DIMENSION(:) :: irecvbuf -real(real64), ALLOCATABLE, DIMENSION(:) :: sendbuf +real(kind=WP), ALLOCATABLE, DIMENSION(:) :: sendbuf integer :: node_size #include "associate_part_def.h" #include "associate_part_ass.h" @@ -1430,7 +1430,7 @@ subroutine broadcast_nod2D(arr2D, arr2Dglobal, partit) sendbuf(i) = arr2Dglobal(irecvbuf(i)) ENDDO - CALL MPI_SEND(sendbuf(1), nTS, MPI_DOUBLE_PRECISION, & + CALL MPI_SEND(sendbuf(1), nTS, MPI_WP, & sender, 2, MPI_COMM_FESOM, MPIerr ) DEALLOCATE(irecvbuf, sendbuf) @@ -1439,7 +1439,7 @@ subroutine broadcast_nod2D(arr2D, arr2Dglobal, partit) CALL MPI_SEND( node_size, 1, MPI_INTEGER, 0, 0, MPI_COMM_FESOM, MPIerr ) CALL MPI_SEND( myList_nod2D(1), node_size, MPI_INTEGER, 0, 1, & MPI_COMM_FESOM, MPIerr ) - CALL MPI_RECV( arr2D(1), node_size, MPI_DOUBLE_PRECISION, 0, & + CALL MPI_RECV( arr2D(1), node_size, MPI_WP, 0, & 2, MPI_COMM_FESOM, status, MPIerr ) ENDIF CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) @@ -1456,12 +1456,12 @@ subroutine broadcast_elem3D(arr3D, arr3Dglobal, partit) type(t_partit), intent(in), target :: partit INTEGER :: nz, counter,nl1 integer :: i, n, nTS, sender, status(MPI_STATUS_SIZE) -real(real64) :: arr3D(:,:) -real(real64) :: arr3Dglobal(:,:) +real(kind=WP) :: arr3D(:,:) +real(kind=WP) :: arr3Dglobal(:,:) integer :: elem_size INTEGER, ALLOCATABLE, DIMENSION(:) :: irecvbuf -real(real64), ALLOCATABLE, DIMENSION(:) :: sendbuf, recvbuf +real(kind=WP), ALLOCATABLE, DIMENSION(:) :: sendbuf, recvbuf #include "associate_part_def.h" #include "associate_part_ass.h" @@ -1488,7 +1488,7 @@ subroutine broadcast_elem3D(arr3D, arr3Dglobal, partit) ENDDO ENDDO - CALL MPI_SEND(sendbuf(1), nTS*nl1, MPI_DOUBLE_PRECISION, & + CALL MPI_SEND(sendbuf(1), nTS*nl1, MPI_WP, & sender, 2, MPI_COMM_FESOM, MPIerr ) DEALLOCATE(irecvbuf, sendbuf) @@ -1499,7 +1499,7 @@ subroutine broadcast_elem3D(arr3D, arr3Dglobal, partit) MPI_COMM_FESOM, MPIerr ) ALLOCATE(recvbuf(elem_size*nl1)) - CALL MPI_RECV( recvbuf(1), elem_size*nl1, MPI_DOUBLE_PRECISION, 0, & + CALL MPI_RECV( recvbuf(1), elem_size*nl1, MPI_WP, 0, & 2, MPI_COMM_FESOM, status, MPIerr ) counter=0 DO n = 1, elem_size @@ -1524,11 +1524,11 @@ subroutine broadcast_elem2D(arr2D, arr2Dglobal, partit) IMPLICIT NONE type(t_partit), intent(in), target :: partit integer :: i, n, nTS, sender, status(MPI_STATUS_SIZE) -real(real64) :: arr2D(:) -real(real64) :: arr2Dglobal(:) +real(kind=WP) :: arr2D(:) +real(kind=WP) :: arr2Dglobal(:) integer :: elem_size INTEGER, ALLOCATABLE, DIMENSION(:) :: irecvbuf -real(real64), ALLOCATABLE, DIMENSION(:) :: sendbuf +real(kind=WP), ALLOCATABLE, DIMENSION(:) :: sendbuf #include "associate_part_def.h" #include "associate_part_ass.h" @@ -1550,7 +1550,7 @@ subroutine broadcast_elem2D(arr2D, arr2Dglobal, partit) sendbuf(i) = arr2Dglobal(irecvbuf(i)) ENDDO - CALL MPI_SEND(sendbuf(1), nTS, MPI_DOUBLE_PRECISION, & + CALL MPI_SEND(sendbuf(1), nTS, MPI_WP, & sender, 2, MPI_COMM_FESOM, MPIerr ) DEALLOCATE(irecvbuf, sendbuf) @@ -1559,7 +1559,7 @@ subroutine broadcast_elem2D(arr2D, arr2Dglobal, partit) CALL MPI_SEND( elem_size, 1, MPI_INTEGER, 0, 0, MPI_COMM_FESOM, MPIerr ) CALL MPI_SEND( myList_elem2D(1), elem_size, MPI_INTEGER, 0, 1, & MPI_COMM_FESOM, MPIerr ) - CALL MPI_RECV( arr2D(1), elem_size, MPI_DOUBLE_PRECISION, 0, & + CALL MPI_RECV( arr2D(1), elem_size, MPI_WP, 0, & 2, MPI_COMM_FESOM, status, MPIerr ) ENDIF CALL MPI_BARRIER(MPI_COMM_FESOM,MPIerr) @@ -2622,11 +2622,11 @@ subroutine gather_edg2D(arr2D, arr2Dglobal, partit) USE MOD_PARSUP IMPLICIT NONE type(t_partit), intent(in), target :: partit -real(real64) :: arr2D(:) -real(real64) :: arr2Dglobal(:) +real(kind=WP) :: arr2D(:) +real(kind=WP) :: arr2Dglobal(:) integer :: i, n, buf_size, sender, status(MPI_STATUS_SIZE) INTEGER, ALLOCATABLE, DIMENSION(:) :: ibuf -REAL(real64), ALLOCATABLE, DIMENSION(:) :: rbuf +real(kind=WP), ALLOCATABLE, DIMENSION(:) :: rbuf #include "associate_part_def.h" #include "associate_part_ass.h" @@ -2640,12 +2640,12 @@ subroutine gather_edg2D(arr2D, arr2Dglobal, partit) !PS CALL MPI_RECV(ibuf(1), buf_size, MPI_INTEGER, sender, & !PS 1, MPI_COMM_FESOM, status, MPIerr ) -!PS CALL MPI_RECV(rbuf(1), buf_size, MPI_DOUBLE_PRECISION, sender, & +!PS CALL MPI_RECV(rbuf(1), buf_size, MPI_WP, sender, & !PS 2, MPI_COMM_FESOM, status, MPIerr ) CALL MPI_RECV(ibuf, buf_size, MPI_INTEGER, sender, & 1, MPI_COMM_FESOM, status, MPIerr ) - CALL MPI_RECV(rbuf, buf_size, MPI_DOUBLE_PRECISION, sender, & + CALL MPI_RECV(rbuf, buf_size, MPI_WP, sender, & 2, MPI_COMM_FESOM, status, MPIerr ) arr2Dglobal(ibuf)=rbuf DEALLOCATE(ibuf, rbuf) @@ -2654,13 +2654,13 @@ subroutine gather_edg2D(arr2D, arr2Dglobal, partit) !PS CALL MPI_SEND( myDim_edge2D, 1, MPI_INTEGER, 0, 0, MPI_COMM_FESOM, MPIerr ) !PS CALL MPI_SEND( myList_edge2D(1), myDim_edge2D, MPI_INTEGER, 0, 1, & !PS MPI_COMM_FESOM, MPIerr ) -!PS CALL MPI_SEND( arr2D(1), myDim_edge2D, MPI_DOUBLE_PRECISION, 0, 2,& +!PS CALL MPI_SEND( arr2D(1), myDim_edge2D, MPI_WP, 0, 2,& !PS MPI_COMM_FESOM, MPIerr ) CALL MPI_SEND( myDim_edge2D, 1, MPI_INTEGER, 0, 0, MPI_COMM_FESOM, MPIerr ) CALL MPI_SEND( myList_edge2D(1:myDim_edge2D), myDim_edge2D, MPI_INTEGER, 0, 1, & MPI_COMM_FESOM, MPIerr ) - CALL MPI_SEND( arr2D(1:myDim_edge2D), myDim_edge2D, MPI_DOUBLE_PRECISION, 0, 2,& + CALL MPI_SEND( arr2D(1:myDim_edge2D), myDim_edge2D, MPI_WP, 0, 2,& MPI_COMM_FESOM, MPIerr ) ENDIF diff --git a/src/gen_ic3d.F90 b/src/gen_ic3d.F90 index 445928bed..837153427 100644 --- a/src/gen_ic3d.F90 +++ b/src/gen_ic3d.F90 @@ -21,6 +21,7 @@ MODULE g_ic3d USE g_support USE g_config, only: dummy, ClimateDataPath, use_cavity USE g_clock, only: r_restart + USE io_netcdf_nf_interface, only: nf_get_vara_x IMPLICIT NONE @@ -192,14 +193,14 @@ SUBROUTINE nc_readGrid(partit) if (partit%mype==0) then nf_start(1)=1 nf_edges(1)=nc_Nlat - iost = nf_get_vara_double(ncid, id_lat, nf_start, nf_edges, nc_lat) + iost = nf_get_vara_x(ncid, id_lat, nf_start, nf_edges, nc_lat) end if call MPI_BCast(iost, 1, MPI_INTEGER, 0, partit%MPI_COMM_FESOM, ierror) call check_nferr(iost,filename,partit) if (partit%mype==0) then nf_start(1)=1 nf_edges(1)=nc_Nlon-2 - iost = nf_get_vara_double(ncid, id_lon, nf_start, nf_edges, nc_lon(2:nc_Nlon-1)) + iost = nf_get_vara_x(ncid, id_lon, nf_start, nf_edges, nc_lon(2:nc_Nlon-1)) nc_lon(1) =nc_lon(nc_Nlon-1) nc_lon(nc_Nlon) =nc_lon(2) end if @@ -209,15 +210,15 @@ SUBROUTINE nc_readGrid(partit) if (partit%mype==0) then nf_start(1)=1 nf_edges(1)=nc_Ndepth - iost = nf_get_vara_double(ncid, id_depth, nf_start, nf_edges,nc_depth) + iost = nf_get_vara_x(ncid, id_depth, nf_start, nf_edges,nc_depth) if (nc_depth(2) < 0.) nc_depth=-nc_depth end if call MPI_BCast(iost, 1, MPI_INTEGER, 0, partit%MPI_COMM_FESOM, ierror) call check_nferr(iost,filename,partit) - call MPI_BCast(nc_lon, nc_Nlon, MPI_DOUBLE_PRECISION, 0, partit%MPI_COMM_FESOM, ierror) - call MPI_BCast(nc_lat, nc_Nlat, MPI_DOUBLE_PRECISION, 0, partit%MPI_COMM_FESOM, ierror) - call MPI_BCast(nc_depth, nc_Ndepth, MPI_DOUBLE_PRECISION, 0, partit%MPI_COMM_FESOM, ierror) + call MPI_BCast(nc_lon, nc_Nlon, MPI_WP, 0, partit%MPI_COMM_FESOM, ierror) + call MPI_BCast(nc_lat, nc_Nlat, MPI_WP, 0, partit%MPI_COMM_FESOM, ierror) + call MPI_BCast(nc_depth, nc_Ndepth, MPI_WP, 0, partit%MPI_COMM_FESOM, ierror) if (partit%mype==0) then iost = nf_close(ncid) @@ -341,11 +342,13 @@ SUBROUTINE getcoeffld(tracers, partit, mesh) integer :: nl1, ul1 real(wp) :: denom, x1, x2, y1, y2, x, y, d1,d2, aux_z real(wp), allocatable, dimension(:,:,:) :: ncdata + real(wp), allocatable, dimension(:) :: ncdata_inner ! rank-1 read buffer (generic nf_get_vara_x only resolves rank-1 actuals) real(wp), allocatable, dimension(:) :: data1d integer :: elnodes(3) integer :: ierror ! return error code integer :: NO_FILL ! 0=no fillval, 1=fillval real(wp) :: FILL_VALUE + real(8) :: FILL_VALUE_r8 ! double temp: file stores fill/missing as double #include "associate_part_def.h" #include "associate_mesh_def.h" #include "associate_part_ass.h" @@ -364,22 +367,24 @@ SUBROUTINE getcoeffld(tracers, partit, mesh) ! get variable id if (mype==0) then iost = nf_inq_varid(ncid, varname, id_data) - iost = nf_inq_var_fill(ncid, id_data, NO_FILL, FILL_VALUE) ! FillValue defined? + iost = nf_inq_var_fill(ncid, id_data, NO_FILL, FILL_VALUE_r8) ! FillValue defined? if (NO_FILL==1) then ! No _FillValue attribute found, try missing_value attribute - iost = nf_get_att_double(ncid, id_data, 'missing_value', FILL_VALUE) + iost = nf_get_att_double(ncid, id_data, 'missing_value', FILL_VALUE_r8) if (iost /= NF_NOERR) then ! Neither _FillValue nor missing_value found, use NetCDF default fill value - FILL_VALUE = NF_FILL_DOUBLE ! 9.9692099683868690e+36 - print *, 'No _FillValue or missing_value in ', trim(filename), ', using NetCDF default:', FILL_VALUE + FILL_VALUE_r8 = NF_FILL_DOUBLE ! 9.9692099683868690e+36 + print *, 'No _FillValue or missing_value in ', trim(filename), ', using NetCDF default:', FILL_VALUE_r8 else - print *, 'Using missing_value from ', trim(filename), ':', FILL_VALUE + print *, 'Using missing_value from ', trim(filename), ':', FILL_VALUE_r8 end if else - print *, 'Using _FillValue from ', trim(filename), ':', FILL_VALUE + print *, 'Using _FillValue from ', trim(filename), ':', FILL_VALUE_r8 end if + ! cast to working precision; comparison below uses the same WP rounding as the data read + FILL_VALUE = real(FILL_VALUE_r8, WP) end if - call MPI_BCast(FILL_VALUE, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(FILL_VALUE, 1, MPI_WP, 0, MPI_COMM_FESOM, ierror) call MPI_BCast(iost, 1, MPI_INTEGER, 0, MPI_COMM_FESOM, ierror) call check_nferr(iost,filename,partit) !read data from file @@ -390,7 +395,12 @@ SUBROUTINE getcoeffld(tracers, partit, mesh) nf_edges(2)=nc_Nlat nf_start(3)=1 nf_edges(3)=nc_Ndepth - iost = nf_get_vara_double(ncid, id_data, nf_start, nf_edges, ncdata(2:nc_Nlon-1,:,:)) + ! read into a contiguous rank-1 WP buffer, then reshape into the strided destination + ! (generic nf_get_vara_x only resolves rank-1 actuals; NetCDF Fortran order has the first dim fastest) + allocate(ncdata_inner((nc_Nlon-2)*nc_Nlat*nc_Ndepth)) + iost = nf_get_vara_x(ncid, id_data, nf_start, nf_edges, ncdata_inner) + ncdata(2:nc_Nlon-1,:,:) = reshape(ncdata_inner, [nc_Nlon-2, nc_Nlat, nc_Ndepth]) + deallocate(ncdata_inner) ncdata(1,:,:) =ncdata(nc_Nlon-1,:,:) ncdata(nc_Nlon,:,:)=ncdata(2,:,:) @@ -410,7 +420,7 @@ SUBROUTINE getcoeffld(tracers, partit, mesh) end if call MPI_BCast(iost, 1, MPI_INTEGER, 0, MPI_COMM_FESOM, ierror) call check_nferr(iost,filename,partit) - call MPI_BCast(ncdata, nc_Nlon*nc_Nlat*nc_Ndepth, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(ncdata, nc_Nlon*nc_Nlat*nc_Ndepth, MPI_WP, 0, MPI_COMM_FESOM, ierror) ! bilinear space interpolation, ! data is assumed to be sampled on a regular grid do ii = 1, myDim_nod2d @@ -628,41 +638,41 @@ SUBROUTINE do_ic3d(tracers, partit, mesh) locO2min = min(locO2min,minval(tracers%data(24)%values(mesh%ulevels_nod2D(n):mesh%nlevels_nod2D(n)-1,n)) ) #endif end do - call MPI_AllREDUCE(locTmax , glo , 1, MPI_DOUBLE_PRECISION, MPI_MAX, partit%MPI_COMM_FESOM, partit%MPIerr) + call MPI_AllREDUCE(locTmax , glo , 1, MPI_WP, MPI_MAX, partit%MPI_COMM_FESOM, partit%MPIerr) if (partit%mype==0) write(*,*) ' |-> gobal max init. temp. =', glo - call MPI_AllREDUCE(locTmin , glo , 1, MPI_DOUBLE_PRECISION, MPI_MIN, partit%MPI_COMM_FESOM, partit%MPIerr) + call MPI_AllREDUCE(locTmin , glo , 1, MPI_WP, MPI_MIN, partit%MPI_COMM_FESOM, partit%MPIerr) if (partit%mype==0) write(*,*) ' |-> gobal min init. temp. =', glo - call MPI_AllREDUCE(locSmax , glo , 1, MPI_DOUBLE_PRECISION, MPI_MAX, partit%MPI_COMM_FESOM, partit%MPIerr) + call MPI_AllREDUCE(locSmax , glo , 1, MPI_WP, MPI_MAX, partit%MPI_COMM_FESOM, partit%MPIerr) if (partit%mype==0) write(*,*) ' |-> gobal max init. salt. =', glo - call MPI_AllREDUCE(locSmin , glo , 1, MPI_DOUBLE_PRECISION, MPI_MIN, partit%MPI_COMM_FESOM, partit%MPIerr) + call MPI_AllREDUCE(locSmin , glo , 1, MPI_WP, MPI_MIN, partit%MPI_COMM_FESOM, partit%MPIerr) if (partit%mype==0) write(*,*) ' `-> gobal min init. salt. =', glo #if defined(__recom) if (partit%mype==0) write(*,*) "Sanity check for REcoM variables" - call MPI_AllREDUCE(locDINmax , glo , 1, MPI_DOUBLE_PRECISION, MPI_MAX, partit%MPI_COMM_FESOM, partit%MPIerr) + call MPI_AllREDUCE(locDINmax , glo , 1, MPI_WP, MPI_MAX, partit%MPI_COMM_FESOM, partit%MPIerr) if (partit%mype==0) write(*,*) ' |-> gobal max init. DIN. =', glo - call MPI_AllREDUCE(locDINmin , glo , 1, MPI_DOUBLE_PRECISION, MPI_MIN, partit%MPI_COMM_FESOM, partit%MPIerr) + call MPI_AllREDUCE(locDINmin , glo , 1, MPI_WP, MPI_MIN, partit%MPI_COMM_FESOM, partit%MPIerr) if (partit%mype==0) write(*,*) ' |-> gobal min init. DIN. =', glo - call MPI_AllREDUCE(locDICmax , glo , 1, MPI_DOUBLE_PRECISION, MPI_MAX, partit%MPI_COMM_FESOM, partit%MPIerr) + call MPI_AllREDUCE(locDICmax , glo , 1, MPI_WP, MPI_MAX, partit%MPI_COMM_FESOM, partit%MPIerr) if (partit%mype==0) write(*,*) ' |-> gobal max init. DIC. =', glo - call MPI_AllREDUCE(locDICmin , glo , 1, MPI_DOUBLE_PRECISION, MPI_MIN, partit%MPI_COMM_FESOM, partit%MPIerr) + call MPI_AllREDUCE(locDICmin , glo , 1, MPI_WP, MPI_MIN, partit%MPI_COMM_FESOM, partit%MPIerr) if (partit%mype==0) write(*,*) ' |-> gobal min init. DIC. =', glo - call MPI_AllREDUCE(locAlkmax , glo , 1, MPI_DOUBLE_PRECISION, MPI_MAX, partit%MPI_COMM_FESOM, partit%MPIerr) + call MPI_AllREDUCE(locAlkmax , glo , 1, MPI_WP, MPI_MAX, partit%MPI_COMM_FESOM, partit%MPIerr) if (partit%mype==0) write(*,*) ' |-> gobal max init. Alk. =', glo - call MPI_AllREDUCE(locAlkmin , glo , 1, MPI_DOUBLE_PRECISION, MPI_MIN, partit%MPI_COMM_FESOM, partit%MPIerr) + call MPI_AllREDUCE(locAlkmin , glo , 1, MPI_WP, MPI_MIN, partit%MPI_COMM_FESOM, partit%MPIerr) if (partit%mype==0) write(*,*) ' |-> gobal min init. Alk. =', glo - call MPI_AllREDUCE(locDSimax , glo , 1, MPI_DOUBLE_PRECISION, MPI_MAX, partit%MPI_COMM_FESOM, partit%MPIerr) + call MPI_AllREDUCE(locDSimax , glo , 1, MPI_WP, MPI_MAX, partit%MPI_COMM_FESOM, partit%MPIerr) if (partit%mype==0) write(*,*) ' |-> gobal max init. DSi. =', glo - call MPI_AllREDUCE(locDSimin , glo , 1, MPI_DOUBLE_PRECISION, MPI_MIN, partit%MPI_COMM_FESOM, partit%MPIerr) + call MPI_AllREDUCE(locDSimin , glo , 1, MPI_WP, MPI_MIN, partit%MPI_COMM_FESOM, partit%MPIerr) if (partit%mype==0) write(*,*) ' |-> gobal min init. DSi. =', glo - call MPI_AllREDUCE(locDFemax , glo , 1, MPI_DOUBLE_PRECISION, MPI_MAX, partit%MPI_COMM_FESOM, partit%MPIerr) + call MPI_AllREDUCE(locDFemax , glo , 1, MPI_WP, MPI_MAX, partit%MPI_COMM_FESOM, partit%MPIerr) if (partit%mype==0) write(*,*) ' |-> gobal max init. DFe. =', glo - call MPI_AllREDUCE(locDFemin , glo , 1, MPI_DOUBLE_PRECISION, MPI_MIN, partit%MPI_COMM_FESOM, partit%MPIerr) + call MPI_AllREDUCE(locDFemin , glo , 1, MPI_WP, MPI_MIN, partit%MPI_COMM_FESOM, partit%MPIerr) if (partit%mype==0) write(*,*) ' `-> gobal min init. DFe. =', glo - call MPI_AllREDUCE(locO2max , glo , 1, MPI_DOUBLE_PRECISION, MPI_MAX, partit%MPI_COMM_FESOM, partit%MPIerr) + call MPI_AllREDUCE(locO2max , glo , 1, MPI_WP, MPI_MAX, partit%MPI_COMM_FESOM, partit%MPIerr) if (partit%mype==0) write(*,*) ' |-> gobal max init. O2. =', glo - call MPI_AllREDUCE(locO2min , glo , 1, MPI_DOUBLE_PRECISION, MPI_MIN, partit%MPI_COMM_FESOM, partit%MPIerr) + call MPI_AllREDUCE(locO2min , glo , 1, MPI_WP, MPI_MIN, partit%MPI_COMM_FESOM, partit%MPIerr) if (partit%mype==0) write(*,*) ' `-> gobal min init. O2. =', glo #endif diff --git a/src/gen_modules_partitioning.F90 b/src/gen_modules_partitioning.F90 index cac1ba210..ad51fffcd 100644 --- a/src/gen_modules_partitioning.F90 +++ b/src/gen_modules_partitioning.F90 @@ -260,7 +260,7 @@ subroutine init_mpi_types(partit, mesh) blocklen_tmp(1:nb) = blocklen(1:nb)*n_val displace_tmp(1:nb) = displace(1:nb)*n_val - call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_DOUBLE_PRECISION, & + call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_WP, & partit%r_mpitype_elem2D(n,n_val), partit%MPIerr) call MPI_TYPE_COMMIT(partit%r_mpitype_elem2D(n,n_val), partit%MPIerr) @@ -276,7 +276,7 @@ subroutine init_mpi_types(partit, mesh) print *,"out of bounds error, lbound:",lbound(partit%r_mpitype_elem3D), "indices:", n,nl1,n_val, "ubound:",ubound(partit%r_mpitype_elem3D), __FILE__,__LINE__ stop 1 end if - call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_DOUBLE_PRECISION, & + call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_WP, & partit%r_mpitype_elem3D(n,nl1,n_val), partit%MPIerr) call MPI_TYPE_COMMIT(partit%r_mpitype_elem3D(n,nl1,n_val), partit%MPIerr) @@ -306,7 +306,7 @@ subroutine init_mpi_types(partit, mesh) blocklen_tmp(1:nb) = blocklen(1:nb)*n_val displace_tmp(1:nb) = displace(1:nb)*n_val - call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_DOUBLE_PRECISION, & + call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_WP, & partit%s_mpitype_elem2D(n, n_val), partit%MPIerr) call MPI_TYPE_COMMIT(partit%s_mpitype_elem2D(n, n_val), partit%MPIerr) @@ -316,7 +316,7 @@ subroutine init_mpi_types(partit, mesh) blocklen_tmp(1:nb) = blocklen(1:nb)*n_val*nl1 displace_tmp(1:nb) = displace(1:nb)*n_val*nl1 - call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_DOUBLE_PRECISION, & + call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_WP, & partit%s_mpitype_elem3D(n,nl1,n_val), partit%MPIerr) call MPI_TYPE_COMMIT(partit%s_mpitype_elem3D(n,nl1,n_val), partit%MPIerr) @@ -347,7 +347,7 @@ subroutine init_mpi_types(partit, mesh) DO n_val=1,4 - call MPI_TYPE_INDEXED(nb, blocklen, displace, MPI_DOUBLE_PRECISION, & + call MPI_TYPE_INDEXED(nb, blocklen, displace, MPI_WP, & partit%r_mpitype_elem2D_full(n,n_val), partit%MPIerr) call MPI_TYPE_COMMIT(partit%r_mpitype_elem2D_full(n, n_val), partit%MPIerr) @@ -362,7 +362,7 @@ subroutine init_mpi_types(partit, mesh) print *,"out of bounds error, lbound:",lbound(partit%r_mpitype_elem3D_full), "indices:", n,nl1,n_val, "ubound:", ubound(partit%r_mpitype_elem3D_full), __FILE__,__LINE__ stop 1 end if - call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_DOUBLE_PRECISION, & + call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_WP, & partit%r_mpitype_elem3D_full(n,nl1,n_val), partit%MPIerr) call MPI_TYPE_COMMIT(partit%r_mpitype_elem3D_full(n,nl1,n_val), partit%MPIerr) @@ -392,7 +392,7 @@ subroutine init_mpi_types(partit, mesh) call MPI_TYPE_COMMIT(partit%s_mpitype_elem2D_full_i(n), partit%MPIerr) DO n_val=1,4 - call MPI_TYPE_INDEXED(nb, blocklen, displace, MPI_DOUBLE_PRECISION, & + call MPI_TYPE_INDEXED(nb, blocklen, displace, MPI_WP, & partit%s_mpitype_elem2D_full(n,n_val), partit%MPIerr) call MPI_TYPE_COMMIT(partit%s_mpitype_elem2D_full(n,n_val), partit%MPIerr) @@ -401,7 +401,7 @@ subroutine init_mpi_types(partit, mesh) blocklen_tmp(1:nb) = blocklen(1:nb)*n_val*nl1 displace_tmp(1:nb) = displace(1:nb)*n_val*nl1 - call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_DOUBLE_PRECISION, & + call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_WP, & partit%s_mpitype_elem3D_full(n,nl1,n_val), partit%MPIerr) call MPI_TYPE_COMMIT(partit%s_mpitype_elem3D_full(n,nl1,n_val), partit%MPIerr) @@ -450,7 +450,7 @@ subroutine init_mpi_types(partit, mesh) endif enddo - call MPI_TYPE_INDEXED(nb, blocklen, displace, MPI_DOUBLE_PRECISION, & + call MPI_TYPE_INDEXED(nb, blocklen, displace, MPI_WP, & partit%r_mpitype_nod2D(n), partit%MPIerr) call MPI_TYPE_INDEXED(nb, blocklen, displace, MPI_INTEGER, & @@ -465,7 +465,7 @@ subroutine init_mpi_types(partit, mesh) blocklen_tmp(1:nb) = blocklen(1:nb)*n_val*nl1 displace_tmp(1:nb) = displace(1:nb)*n_val*nl1 - call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_DOUBLE_PRECISION, & + call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_WP, & partit%r_mpitype_nod3D(n,nl1,n_val), partit%MPIerr) call MPI_TYPE_COMMIT(partit%r_mpitype_nod3D(n,nl1,n_val), partit%MPIerr) @@ -490,7 +490,7 @@ subroutine init_mpi_types(partit, mesh) endif enddo - call MPI_TYPE_INDEXED(nb, blocklen, displace, MPI_DOUBLE_PRECISION, & + call MPI_TYPE_INDEXED(nb, blocklen, displace, MPI_WP, & partit%s_mpitype_nod2D(n), partit%MPIerr) call MPI_TYPE_INDEXED(nb, blocklen, displace, MPI_INTEGER, & @@ -505,7 +505,7 @@ subroutine init_mpi_types(partit, mesh) blocklen_tmp(1:nb) = blocklen(1:nb)*n_val*nl1 displace_tmp(1:nb) = displace(1:nb)*n_val*nl1 - call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_DOUBLE_PRECISION, & + call MPI_TYPE_INDEXED(nb, blocklen_tmp, displace_tmp, MPI_WP, & partit%s_mpitype_nod3D(n,nl1,n_val), partit%MPIerr) call MPI_TYPE_COMMIT(partit%s_mpitype_nod3D(n,nl1,n_val), partit%MPIerr) diff --git a/src/gen_modules_read_NetCDF.F90 b/src/gen_modules_read_NetCDF.F90 index 8ab8b94fd..3bc413d28 100755 --- a/src/gen_modules_read_NetCDF.F90 +++ b/src/gen_modules_read_NetCDF.F90 @@ -28,11 +28,11 @@ subroutine read_other_NetCDF(file, vari, itime, model_2Darray, check_dummy, do_o integer :: status, ncid, varid integer :: lonid, latid integer :: istart(3), icount(3), elnodes(3) - real(real64) :: x, y, miss, aux, xmin, elnodes_x(3) - real(real64), allocatable :: lon(:), lat(:) - real(real64), allocatable :: ncdata(:,:), ncdata_temp(:,:) - real(real64), allocatable :: temp_x(:), temp_y(:) - real(real64) :: model_2Darray(partit%myDim_nod2d+partit%eDim_nod2D) + real(kind=WP) :: x, y, miss, aux, xmin, elnodes_x(3) + real(kind=WP), allocatable :: lon(:), lat(:) + real(kind=WP), allocatable :: ncdata(:,:), ncdata_temp(:,:) + real(kind=WP), allocatable :: temp_x(:), temp_y(:) + real(kind=WP) :: model_2Darray(partit%myDim_nod2d+partit%eDim_nod2D) character(*) :: vari character(*) :: file logical :: check_dummy, do_onvert @@ -74,7 +74,7 @@ subroutine read_other_NetCDF(file, vari, itime, model_2Darray, check_dummy, do_o status=nf90_inq_varid(ncid, 'lat', varid) status=nf90_get_var(ncid, varid, lat, start=(/1/), count=(/latlen/)) end if - call MPI_BCast(lat, latlen, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(lat, latlen, MPI_WP, 0, MPI_COMM_FESOM, ierror) ! lon allocate(lon(lonlen)) @@ -82,7 +82,7 @@ subroutine read_other_NetCDF(file, vari, itime, model_2Darray, check_dummy, do_o status=nf90_inq_varid(ncid, 'lon', varid) status=nf90_get_var(ncid, varid, lon, start=(/1/), count=(/lonlen/)) end if - call MPI_BCast(lon, lonlen, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(lon, lonlen, MPI_WP, 0, MPI_COMM_FESOM, ierror) ! make sure range 0. - 360. do n=1,lonlen @@ -106,8 +106,8 @@ subroutine read_other_NetCDF(file, vari, itime, model_2Darray, check_dummy, do_o ! close file status=nf90_close(ncid) end if - call MPI_BCast(ncdata, lonlen*latlen, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) - call MPI_BCast(miss, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(ncdata, lonlen*latlen, MPI_WP, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(miss, 1, MPI_WP, 0, MPI_COMM_FESOM, ierror) !write(*,*)'miss', miss !write(*,*)'raw',minval(ncdata),maxval(ncdata) ncdata_temp=ncdata @@ -209,11 +209,11 @@ subroutine read_surf_hydrography_NetCDF(file, vari, itime, model_2Darray, partit integer :: status, ncid, varid integer :: lonid, latid, drain_num integer :: istart(4), icount(4) - real(real64) :: x, y, miss - real(real64), allocatable :: lon(:), lat(:) - real(real64), allocatable :: ncdata(:,:) - real(real64), allocatable :: temp_x(:), temp_y(:) - real(real64) :: model_2Darray(partit%myDim_nod2d+partit%eDim_nod2D) + real(kind=WP) :: x, y, miss + real(kind=WP), allocatable :: lon(:), lat(:) + real(kind=WP), allocatable :: ncdata(:,:) + real(kind=WP), allocatable :: temp_x(:), temp_y(:) + real(kind=WP) :: model_2Darray(partit%myDim_nod2d+partit%eDim_nod2D) character(15) :: vari character(300) :: file logical :: check_dummy @@ -252,7 +252,7 @@ subroutine read_surf_hydrography_NetCDF(file, vari, itime, model_2Darray, partit status=nf90_inq_varid(ncid, 'lat', varid) status=nf90_get_var(ncid, varid, lat, start=(/1/), count=(/latlen/)) end if - call MPI_BCast(lat, latlen, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(lat, latlen, MPI_WP, 0, MPI_COMM_FESOM, ierror) ! lon allocate(lon(lonlen)) @@ -260,7 +260,7 @@ subroutine read_surf_hydrography_NetCDF(file, vari, itime, model_2Darray, partit status=nf90_inq_varid(ncid, 'lon', varid) status=nf90_get_var(ncid, varid, lon, start=(/1/), count=(/lonlen/)) end if - call MPI_BCast(lon, lonlen, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(lon, lonlen, MPI_WP, 0, MPI_COMM_FESOM, ierror) ! make sure range 0. - 360. do n=1,lonlen @@ -286,8 +286,8 @@ subroutine read_surf_hydrography_NetCDF(file, vari, itime, model_2Darray, partit !close file status=nf90_close(ncid) end if - call MPI_BCast(ncdata, lonlen*latlen, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) - call MPI_BCast(miss, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(ncdata, lonlen*latlen, MPI_WP, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(miss, 1, MPI_WP, 0, MPI_COMM_FESOM, ierror) ! the next step is to interpolate data to model grids ! model grid coordinates num=myDim_nod2d+eDim_nod2d @@ -331,8 +331,8 @@ subroutine read_2ddata_on_grid_NetCDF(file, vari, itime, model_2Darray, partit, integer :: itime integer :: status, ncid, varid integer :: istart(2), icount(2) - real(real64) :: ncdata(mesh%nod2D) - real(real64), intent(out) :: model_2Darray(partit%myDim_nod2D+partit%eDim_nod2D) + real(kind=WP) :: ncdata(mesh%nod2D) + real(kind=WP), intent(out) :: model_2Darray(partit%myDim_nod2D+partit%eDim_nod2D) character(*), intent(in) :: file character(*), intent(in) :: vari integer :: ierror ! return error code @@ -362,7 +362,7 @@ subroutine read_2ddata_on_grid_NetCDF(file, vari, itime, model_2Darray, partit, status=nf90_get_var(ncid, varid, ncdata, start=istart, count=icount) status=nf90_close(ncid) end if - call MPI_BCast(ncdata, nod2D, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(ncdata, nod2D, MPI_WP, 0, MPI_COMM_FESOM, ierror) model_2Darray=ncdata(myList_nod2D) end subroutine read_2ddata_on_grid_NetCDF end module g_read_other_NetCDF diff --git a/src/gen_support.F90 b/src/gen_support.F90 index 01ef9c0c6..960987296 100644 --- a/src/gen_support.F90 +++ b/src/gen_support.F90 @@ -346,7 +346,7 @@ subroutine integrate_nod_2D(data, int2D, partit, mesh) !$OMP END PARALLEL #endif int2D=0.0_WP - call MPI_AllREDUCE(lval, int2D, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & + call MPI_AllREDUCE(lval, int2D, 1, MPI_WP, MPI_SUM, & MPI_COMM_FESOM, MPIerr) end subroutine integrate_nod_2D ! @@ -391,7 +391,7 @@ subroutine integrate_nod_3D(data, int3D, partit, mesh) !$OMP END PARALLEL DO int3D=0.0_WP - call MPI_AllREDUCE(lval, int3D, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & + call MPI_AllREDUCE(lval, int3D, 1, MPI_WP, MPI_SUM, & MPI_COMM_FESOM, MPIerr) end subroutine integrate_nod_3D ! @@ -421,7 +421,7 @@ subroutine extrap_nod3D(arr, partit, mesh) !___________________________________________________________________________ loc_max=maxval(arr(1,:)) glob_max=0._WP - call MPI_AllREDUCE(loc_max, glob_max, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc_max, glob_max, 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) glob_sum=-1 !___________________________________________________________________________ @@ -485,7 +485,7 @@ subroutine extrap_nod3D(arr, partit, mesh) !_______________________________________________________________________ loc_max=maxval(arr(1,:)) - call MPI_AllREDUCE(loc_max, glob_max, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc_max, glob_max, 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) END DO ! --> DO WHILE (glob_max>0.99_WP*dummy) @@ -681,7 +681,7 @@ subroutine integrate_elem_3D(data, int3D, partit, mesh) !$OMP END PARALLEL DO int3D=0.0_WP - call MPI_AllREDUCE(lval, int3D, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & + call MPI_AllREDUCE(lval, int3D, 1, MPI_WP, MPI_SUM, & MPI_COMM_FESOM, MPIerr) end subroutine integrate_elem_3D ! @@ -721,7 +721,7 @@ subroutine integrate_elem_2D(data, int2D, partit, mesh) !$OMP END PARALLEL DO int2D=0.0_WP - call MPI_AllREDUCE(lval, int2D, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & + call MPI_AllREDUCE(lval, int2D, 1, MPI_WP, MPI_SUM, & MPI_COMM_FESOM, MPIerr) end subroutine integrate_elem_2D diff --git a/src/gen_surface_forcing.F90 b/src/gen_surface_forcing.F90 index efb4248bc..894b73559 100644 --- a/src/gen_surface_forcing.F90 +++ b/src/gen_surface_forcing.F90 @@ -387,7 +387,7 @@ SUBROUTINE nc_readTimeGrid(flf, partit) iost = nf90_get_var(ncid, id_time, flf%nc_time, start=(/1/), count=(/flf%nc_Ntime/)) ! digg for calendar attribute in time axis variable end if - call MPI_BCast(flf%nc_time, flf%nc_Ntime, MPI_DOUBLE_PRECISION, 0, partit%MPI_COMM_FESOM, ierror) + call MPI_BCast(flf%nc_time, flf%nc_Ntime, MPI_WP, 0, partit%MPI_COMM_FESOM, ierror) call MPI_BCast(iost, 1, MPI_INTEGER, 0, partit%MPI_COMM_FESOM, ierror) call check_nferr(iost,flf%file_name,partit) @@ -511,8 +511,8 @@ SUBROUTINE nc_readTimeGrid(flf, partit) flf%nc_time(flf%nc_Ntime) = flf%nc_time(flf%nc_Ntime) + (flf%nc_time(flf%nc_Ntime) - flf%nc_time(flf%nc_Ntime-1))/2.0 end if end if - call MPI_BCast(flf%nc_lon, flf%nc_Nlon, MPI_DOUBLE_PRECISION, 0, partit%MPI_COMM_FESOM, ierror) - call MPI_BCast(flf%nc_lat, flf%nc_Nlat, MPI_DOUBLE_PRECISION, 0, partit%MPI_COMM_FESOM, ierror) + call MPI_BCast(flf%nc_lon, flf%nc_Nlon, MPI_WP, 0, partit%MPI_COMM_FESOM, ierror) + call MPI_BCast(flf%nc_lat, flf%nc_Nlat, MPI_WP, 0, partit%MPI_COMM_FESOM, ierror) !___________________________________________________________________________ !flip lat and data in case of lat from -90 to 90 @@ -3229,8 +3229,8 @@ subroutine read_runoff_mapper(file, vari, R, partit, mesh) end do deallocate(ncdata, lon, lat) end if - call MPI_BCast(lon_sparse, number_arrival_points, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) - call MPI_BCast(lat_sparse, number_arrival_points, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(lon_sparse, number_arrival_points, MPI_WP, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(lat_sparse, number_arrival_points, MPI_WP, 0, MPI_COMM_FESOM, ierror) call MPI_BCast(data_sparse, number_arrival_points, MPI_INTEGER, 0, MPI_COMM_FESOM, ierror) drain_num=maxval(data_sparse) ALLOCATE(arrival_area(drain_num)) @@ -3255,7 +3255,7 @@ subroutine read_runoff_mapper(file, vari, R, partit, mesh) do i=1, number_arrival_points dist_min_glo(i)=dist_min(i) end do - call MPI_AllREDUCE(MPI_IN_PLACE , dist_min_glo , number_arrival_points, MPI_DOUBLE, MPI_MIN, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(MPI_IN_PLACE , dist_min_glo , number_arrival_points, MPI_WP, MPI_MIN, MPI_COMM_FESOM, MPIerr) lon_sparse=0.0_WP lat_sparse=0.0_WP @@ -3268,8 +3268,8 @@ subroutine read_runoff_mapper(file, vari, R, partit, mesh) status=status+1 end if end do - call MPI_AllREDUCE(MPI_IN_PLACE , lon_sparse , number_arrival_points, MPI_DOUBLE, MPI_SUM, MPI_COMM_FESOM, MPIerr) - call MPI_AllREDUCE(MPI_IN_PLACE , lat_sparse , number_arrival_points, MPI_DOUBLE, MPI_SUM, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(MPI_IN_PLACE , lon_sparse , number_arrival_points, MPI_WP, MPI_SUM, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(MPI_IN_PLACE , lat_sparse , number_arrival_points, MPI_WP, MPI_SUM, MPI_COMM_FESOM, MPIerr) call MPI_AllREDUCE(MPI_IN_PLACE , status , 1, MPI_INTEGER, MPI_SUM, MPI_COMM_FESOM, MPIerr) if (status/=number_arrival_points) then @@ -3326,7 +3326,7 @@ subroutine read_runoff_mapper(file, vari, R, partit, mesh) end if END DO - call MPI_AllREDUCE(MPI_IN_PLACE , arrival_area, drain_num, MPI_DOUBLE, MPI_SUM, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(MPI_IN_PLACE , arrival_area, drain_num, MPI_WP, MPI_SUM, MPI_COMM_FESOM, MPIerr) DO i=1, drain_num where (RUNOFF_MAPPER%colind==i) diff --git a/src/icb_allocate.F90 b/src/icb_allocate.F90 index 6f0c437a1..5663d5f00 100644 --- a/src/icb_allocate.F90 +++ b/src/icb_allocate.F90 @@ -135,7 +135,7 @@ subroutine allocate_icb(partit, mesh) allocate(elem_area_glob(elem2D)) elem_area_glob=0.0 call gather_elem(elem_area(1:myDim_elem2D), elem_area_glob, partit) - call MPI_Bcast(elem_area_glob, elem2D, MPI_DOUBLE, 0, MPI_COMM_FESOM, MPIERR) + call MPI_Bcast(elem_area_glob, elem2D, MPI_WP, 0, MPI_COMM_FESOM, MPIERR) allocate(vl_block(4*ib_num)) allocate(buoy_props(ib_num,14)) diff --git a/src/icb_elem.F90 b/src/icb_elem.F90 index c0a0c8421..3f355e641 100644 --- a/src/icb_elem.F90 +++ b/src/icb_elem.F90 @@ -712,15 +712,15 @@ end subroutine com_integer ! CALL MPI_RECV(he_has_element, 1, MPI_LOGICAL, MPI_ANY_SOURCE, 0, MPI_COMM_FESOM, status, MPIerr ) ! sender = status(MPI_SOURCE) ! if (he_has_element) then -! CALL MPI_RECV(arr_r, 15, MPI_DOUBLE_PRECISION, sender, 1, MPI_COMM_FESOM, status, MPIerr ) -! CALL MPI_RECV(iceberg_element, 1, MPI_DOUBLE_PRECISION, sender, 2, MPI_COMM_FESOM, status, MPIerr ) +! CALL MPI_RECV(arr_r, 15, MPI_WP, sender, 1, MPI_COMM_FESOM, status, MPIerr ) +! CALL MPI_RECV(iceberg_element, 1, MPI_WP, sender, 2, MPI_COMM_FESOM, status, MPIerr ) ! arr=arr_r ! end if ! end do ! else ! CALL MPI_SEND(i_have_element, 1, MPI_LOGICAL, 0, 0, MPI_COMM_FESOM, MPIerr ) ! if (i_have_element) then -! CALL MPI_SEND(arr, 15, MPI_DOUBLE_PRECISION,0, 1, MPI_COMM_FESOM, MPIerr ) +! CALL MPI_SEND(arr, 15, MPI_WP,0, 1, MPI_COMM_FESOM, MPIerr ) ! CALL MPI_SEND(iceberg_element, 1, MPI_INTEGER,0, 2, MPI_COMM_FESOM, MPIerr ) ! end if ! end if @@ -730,15 +730,15 @@ end subroutine com_integer !! CALL MPI_RECV(he_has_element, 1, MPI_LOGICAL, MPI_ANY_SOURCE, 0, MPI_COMM_WORLD, status, MPIerr ) !! sender = status(MPI_SOURCE) !! if (he_has_element) then -!! CALL MPI_RECV(arr_r, 15, MPI_DOUBLE_PRECISION, sender, 1, MPI_COMM_WORLD, status, MPIerr ) -!! CALL MPI_RECV(iceberg_element, 1, MPI_DOUBLE_PRECISION, sender, 2, MPI_COMM_WORLD, status, MPIerr ) +!! CALL MPI_RECV(arr_r, 15, MPI_WP, sender, 1, MPI_COMM_WORLD, status, MPIerr ) +!! CALL MPI_RECV(iceberg_element, 1, MPI_WP, sender, 2, MPI_COMM_WORLD, status, MPIerr ) !! arr=arr_r !! end if !! end do !! else !! CALL MPI_SEND(i_have_element, 1, MPI_LOGICAL, 0, 0, MPI_COMM_WORLD, MPIerr ) !! if (i_have_element) then -!! CALL MPI_SEND(arr, 15, MPI_DOUBLE_PRECISION,0, 1, MPI_COMM_WORLD, MPIerr ) +!! CALL MPI_SEND(arr, 15, MPI_WP,0, 1, MPI_COMM_WORLD, MPIerr ) !! CALL MPI_SEND(iceberg_element, 1, MPI_INTEGER,0, 2, MPI_COMM_WORLD, MPIerr ) !! end if !! end if @@ -750,9 +750,9 @@ end subroutine com_integer ! !3. datatype - Datentyp der Pufferelemente (handle) ! !4. root - Wurzelproze�; der, welcher sendet (integer) ! !5. comm - Kommunikator (handle) -! CALL MPI_BCAST(arr, 15, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, MPIerr) +! CALL MPI_BCAST(arr, 15, MPI_WP, 0, MPI_COMM_FESOM, MPIerr) ! CALL MPI_BCAST(iceberg_element, 1, MPI_INTEGER, 0, MPI_COMM_FESOM, MPIerr) -! !CALL MPI_BCAST(arr, 15, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, MPIerr) +! !CALL MPI_BCAST(arr, 15, MPI_WP, 0, MPI_COMM_WORLD, MPIerr) ! !CALL MPI_BCAST(iceberg_element, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, MPIerr) ! ! ! kh 10.02.21 @@ -864,9 +864,9 @@ subroutine matrix_inverse_2x2 (A, AINV, DET, elem, coords) integer :: elem REAL, DIMENSION(2), INTENT(IN) :: coords - real(kind=8), dimension(2,2), intent(IN) :: A - real(kind=8), dimension(2,2), intent(OUT) :: AINV - real(kind=8), intent(OUT) :: DET + real(kind=WP), dimension(2,2), intent(IN) :: A + real(kind=WP), dimension(2,2), intent(OUT) :: AINV + real(kind=WP), intent(OUT) :: DET integer :: i,j DET = A(1,1)*A(2,2) - A(1,2)*A(2,1) diff --git a/src/icb_modules.F90 b/src/icb_modules.F90 index fe2f89c2e..db8bd04ff 100644 --- a/src/icb_modules.F90 +++ b/src/icb_modules.F90 @@ -116,7 +116,7 @@ module iceberg_params real,dimension(:), allocatable:: arr_block integer,dimension(:), allocatable:: elem_block integer,dimension(:), allocatable:: pe_block - real(real64), dimension(:), allocatable:: elem_area_glob + real(kind=WP), dimension(:), allocatable:: elem_area_glob real,dimension(:), allocatable:: vl_block !array for output in netcdf diff --git a/src/icb_step.F90 b/src/icb_step.F90 index b1e4b4e78..d238f8f2f 100644 --- a/src/icb_step.F90 +++ b/src/icb_step.F90 @@ -52,7 +52,7 @@ subroutine iceberg_calculation(ice, mesh, partit, dynamics, istep) integer :: istep_end_synced integer:: req, status(MPI_STATUS_SIZE) logical:: completed - real(kind=8) :: t0, t1, t2, t3, t4, t0_restart, t1_restart != + real(kind=WP) :: t0, t1, t2, t3, t4, t0_restart, t1_restart != logical :: firstcall=.true. != logical :: lastsubstep != @@ -148,7 +148,7 @@ subroutine iceberg_calculation(ice, mesh, partit, dynamics, istep) vl_block_red = 0.0 !$omp critical - call MPI_IAllREDUCE(arr_block, arr_block_red, 16*ib_num, MPI_DOUBLE_PRECISION, MPI_SUM, partit%MPI_COMM_FESOM_IB, req, partit%MPIERR_IB) + call MPI_IAllREDUCE(arr_block, arr_block_red, 16*ib_num, MPI_WP, MPI_SUM, partit%MPI_COMM_FESOM_IB, req, partit%MPIERR_IB) !$omp end critical completed = .false. @@ -193,7 +193,7 @@ subroutine iceberg_calculation(ice, mesh, partit, dynamics, istep) !$omp critical - call MPI_IAllREDUCE(vl_block, vl_block_red, 4*ib_num, MPI_DOUBLE_PRECISION, MPI_SUM, partit%MPI_COMM_FESOM_IB, req, partit%MPIERR_IB) + call MPI_IAllREDUCE(vl_block, vl_block_red, 4*ib_num, MPI_WP, MPI_SUM, partit%MPI_COMM_FESOM_IB, req, partit%MPIERR_IB) !$omp end critical completed = .false. @@ -346,7 +346,7 @@ subroutine iceberg_step1(ice, mesh, partit, dynamics, ib, height_ib_single,lengt logical :: i_have_element != real :: left_mype != integer :: old_element != - real(kind=8) :: t0, t1, t2, t3, t4, t5, t6, t7, t8 != + real(kind=WP) :: t0, t1, t2, t3, t4, t5, t6, t7, t8 != != !for restart != logical, save :: firstcall=.true. != @@ -710,8 +710,8 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib_single integer status(MPI_STATUS_SIZE) integer :: num_ib_in_elem, idx real :: area_ib_tot - !real(real64), dimension(:), allocatable :: rbuffer, local_elem_area - real(real64) :: elem_area_tmp + !real(kind=WP), dimension(:), allocatable :: rbuffer, local_elem_area + real(kind=WP) :: elem_area_tmp !iceberg output character :: ib_char*10 @@ -725,7 +725,7 @@ subroutine iceberg_step2(mesh, partit,arr, elem_from_block, ib, height_ib_single logical :: i_have_element != real :: left_mype != integer :: old_element != - real(kind=8) :: t0, t1, t2, t3, t4 != + real(kind=WP) :: t0, t1, t2, t3, t4 != != !for restart != logical, save :: firstcall=.true. != @@ -1981,7 +1981,7 @@ subroutine write_buoy_props_netcdf(partit) integer :: height_id, length_id, width_id integer :: bvl_id, lvlv_id, lvle_id, lvlb_id, felem_id, grounded_id integer :: start(2), count(2) - real(kind=8) :: sec_in_year + real(kind=WP) :: sec_in_year type(t_partit), intent(inout), target :: partit !type(t_ice), intent(inout), target :: ice #include "associate_part_def.h" diff --git a/src/icb_thermo.F90 b/src/icb_thermo.F90 index b2cf375f4..424cadf89 100644 --- a/src/icb_thermo.F90 +++ b/src/icb_thermo.F90 @@ -53,7 +53,7 @@ subroutine iceberg_meltrates(partit, mesh, M_b, M_v, M_e, M_bv, & implicit none ! LA: include latent heat 2023-04-04 - real(kind=8),parameter :: L = 334000. ! [J/Kg] + real(kind=WP),parameter :: L = 334000. ! [J/Kg] real, intent(IN) :: u_ib,v_ib, uo_ib,vo_ib,ua_ib,va_ib !iceberg velo, (int.) ocean & atm velo real, intent(IN) :: uo_keel_ib, vo_keel_ib !ocean velo at iceberg's draft @@ -220,7 +220,7 @@ subroutine iceberg_newdimensions(partit, ib, depth_ib,height_ib,length_ib,width_ real, dimension(4) :: arr integer :: istep ! LA: include latent heat 2023-04-04 - real(kind=8),parameter :: L = 334000. ! [J/Kg] + real(kind=WP),parameter :: L = 334000. ! [J/Kg] type(t_partit), intent(inout), target :: partit #include "associate_part_def.h" @@ -407,42 +407,42 @@ subroutine iceberg_heat_water_fluxes_3eq(partit, ib, M_b, H_b, T_ib,S_ib,v_rel, implicit none integer, INTENT(IN) :: ib - real(kind=8),INTENT(OUT) :: M_b, H_b, t_freeze - real(kind=8),INTENT(IN) :: T_ib, S_ib ! ocean temperature & salinity (at depth 'depth_ib') - real(kind=8),INTENT(IN) :: v_rel, depth_ib ! relative velocity iceberg-ocean (at depth 'depth_ib') - - real (kind=8) :: temp,sal,tin,zice - real (kind=8) :: rhow, rhor, rho - real (kind=8) :: gats1, gats2, gas, gat - real (kind=8) :: ep1,ep2,ep3,ep4,ep5,ep31 - real (kind=8) :: ex1,ex2,ex3,ex4,ex5,ex6 - real (kind=8) :: vt1,sr1,sr2,sf1,sf2,tf1,tf2,tf,sf,seta,re + real(kind=WP),INTENT(OUT) :: M_b, H_b, t_freeze + real(kind=WP),INTENT(IN) :: T_ib, S_ib ! ocean temperature & salinity (at depth 'depth_ib') + real(kind=WP),INTENT(IN) :: v_rel, depth_ib ! relative velocity iceberg-ocean (at depth 'depth_ib') + + real(kind=WP) :: temp,sal,tin,zice + real(kind=WP) :: rhow, rhor, rho + real(kind=WP) :: gats1, gats2, gas, gat + real(kind=WP) :: ep1,ep2,ep3,ep4,ep5,ep31 + real(kind=WP) :: ex1,ex2,ex3,ex4,ex5,ex6 + real(kind=WP) :: vt1,sr1,sr2,sf1,sf2,tf1,tf2,tf,sf,seta,re integer :: n, n3, nk - real(kind=8),parameter :: rp = 0. !reference pressure - real(kind=8),parameter :: a = -0.0575 !Foldvik&Kvinge (1974) - real(kind=8),parameter :: b = 0.0901 - real(kind=8),parameter :: c = 7.61e-4 - - real(kind=8),parameter :: pr = 13.8 !Prandtl number [dimensionless] - real(kind=8),parameter :: sc = 2432. !Schmidt number [dimensionless] - real(kind=8),parameter :: ak = 2.50e-3 !dimensionless drag coeff. - real(kind=8),parameter :: sak1= sqrt(ak) - real(kind=8),parameter :: un = 1.95e-6 !kinematic viscosity [m2/s] - real(kind=8),parameter :: pr1 = pr**(2./3.) !Jenkins (1991) - real(kind=8),parameter :: sc1 = sc**(2./3.) - - real(kind=8),parameter :: tob= -20. !temperatur at the ice surface - !real(kind=8),parameter :: rhoi= 920. !mean ice density - !real(kind=8),parameter :: rhoh2o= 1027.5 !water density - real(kind=8),parameter :: rhoi= 850.0 !mean ice(berg) density (see values in icb_modules.F90) - real(kind=8),parameter :: cpw = 4180.0 !Barnier et al. (1995) - real(kind=8),parameter :: lhf = 3.33e+5 !latent heat of fusion - real(kind=8),parameter :: tdif= 1.54e-6 !thermal conductivity of ice shelf !RG4190 / RG44027 - real(kind=8),parameter :: atk = 273.15 !0 deg C in Kelvin - real(kind=8),parameter :: cpi = 152.5+7.122*(atk+tob) !Paterson:"The Physics of Glaciers" - - real(kind=8),parameter :: L = 334000. ! [J/Kg] + real(kind=WP),parameter :: rp = 0. !reference pressure + real(kind=WP),parameter :: a = -0.0575 !Foldvik&Kvinge (1974) + real(kind=WP),parameter :: b = 0.0901 + real(kind=WP),parameter :: c = 7.61e-4 + + real(kind=WP),parameter :: pr = 13.8 !Prandtl number [dimensionless] + real(kind=WP),parameter :: sc = 2432. !Schmidt number [dimensionless] + real(kind=WP),parameter :: ak = 2.50e-3 !dimensionless drag coeff. + real(kind=WP),parameter :: sak1= sqrt(ak) + real(kind=WP),parameter :: un = 1.95e-6 !kinematic viscosity [m2/s] + real(kind=WP),parameter :: pr1 = pr**(2./3.) !Jenkins (1991) + real(kind=WP),parameter :: sc1 = sc**(2./3.) + + real(kind=WP),parameter :: tob= -20. !temperatur at the ice surface + !real(kind=WP),parameter :: rhoi= 920. !mean ice density + !real(kind=WP),parameter :: rhoh2o= 1027.5 !water density + real(kind=WP),parameter :: rhoi= 850.0 !mean ice(berg) density (see values in icb_modules.F90) + real(kind=WP),parameter :: cpw = 4180.0 !Barnier et al. (1995) + real(kind=WP),parameter :: lhf = 3.33e+5 !latent heat of fusion + real(kind=WP),parameter :: tdif= 1.54e-6 !thermal conductivity of ice shelf !RG4190 / RG44027 + real(kind=WP),parameter :: atk = 273.15 !0 deg C in Kelvin + real(kind=WP),parameter :: cpi = 152.5+7.122*(atk+tob) !Paterson:"The Physics of Glaciers" + + real(kind=WP),parameter :: L = 334000. ! [J/Kg] type(t_partit), intent(inout), target :: partit !==================== MODULES & DECLARATIONS ==========================!= #include "associate_part_def.h" @@ -631,9 +631,9 @@ subroutine fcn_density(t,s,z,rho) use o_PARAM implicit none - real(kind=8), intent(IN) :: t, s, z - real(kind=8), intent(OUT) :: rho - real(kind=8) :: rhopot, bulk + real(kind=WP), intent(IN) :: t, s, z + real(kind=WP), intent(OUT) :: rho + real(kind=WP) :: rhopot, bulk bulk = 19092.56 + t*(209.8925 & - t*(3.041638 - t*(-1.852732e-3 & diff --git a/src/io_fesom_file.F90 b/src/io_fesom_file.F90 index 1706f0cec..b75a5077f 100644 --- a/src/io_fesom_file.F90 +++ b/src/io_fesom_file.F90 @@ -10,8 +10,8 @@ module io_fesom_file_module type var_info integer var_index - real(kind=8), pointer :: external_local_data_ptr(:,:) => null() - real(kind=8), allocatable, dimension(:,:) :: local_data_copy + real(kind=WP), pointer :: external_local_data_ptr(:,:) => null() + real(kind=WP), allocatable, dimension(:,:) :: local_data_copy real(kind=8), allocatable :: global_level_data(:) integer :: global_level_data_size = 0 logical is_elem_based @@ -485,9 +485,9 @@ subroutine specify_node_var_2d(this, name, longname, units, local_data) class(fesom_file_type), intent(inout) :: this character(len=*), intent(in) :: name character(len=*), intent(in) :: units, longname - real(kind=8), target, intent(inout) :: local_data(:) ! todo: be able to set precision + real(kind=WP), target, intent(inout) :: local_data(:) ! todo: be able to set precision ! EO parameters - real(8), pointer :: external_local_data_ptr(:,:) + real(kind=WP), pointer :: external_local_data_ptr(:,:) type(dim_info) level_diminfo !PS write(*,*) "--> specify_node_var_2d:", __LINE__, __FILE__ @@ -502,7 +502,7 @@ subroutine specify_node_var_2dicepack(this, name, longname, units, local_data, n class(fesom_file_type), intent(inout) :: this character(len=*), intent(in) :: name character(len=*), intent(in) :: units, longname - real(kind=8), target, intent(inout) :: local_data(:,:) + real(kind=WP), target, intent(inout) :: local_data(:,:) integer, intent(in) :: ncat! todo: be able to set precision ! EO parameters type(dim_info) level_diminfo, ncat_diminfo @@ -521,7 +521,7 @@ subroutine specify_node_var_3d(this, name, longname, units, local_data) class(fesom_file_type), intent(inout) :: this character(len=*), intent(in) :: name character(len=*), intent(in) :: units, longname - real(kind=8), target, intent(inout) :: local_data(:,:) ! todo: be able to set precision + real(kind=WP), target, intent(inout) :: local_data(:,:) ! todo: be able to set precision ! EO parameters type(dim_info) level_diminfo, depth_diminfo @@ -538,9 +538,9 @@ subroutine specify_elem_var_2d(this, name, longname, units, local_data) class(fesom_file_type), intent(inout) :: this character(len=*), intent(in) :: name character(len=*), intent(in) :: units, longname - real(kind=8), target, intent(inout) :: local_data(:) ! todo: be able to set precision + real(kind=WP), target, intent(inout) :: local_data(:) ! todo: be able to set precision ! EO parameters - real(8), pointer :: external_local_data_ptr(:,:) + real(kind=WP), pointer :: external_local_data_ptr(:,:) type(dim_info) level_diminfo !PS write(*,*) "--> specify_elem_var_2d:", __LINE__, __FILE__ @@ -556,7 +556,7 @@ subroutine specify_elem_var_3d(this, name, longname, units, local_data) class(fesom_file_type), intent(inout) :: this character(len=*), intent(in) :: name character(len=*), intent(in) :: units, longname - real(kind=8), target, intent(inout) :: local_data(:,:) ! todo: be able to set precision + real(kind=WP), target, intent(inout) :: local_data(:,:) ! todo: be able to set precision ! EO parameters type(dim_info) level_diminfo, depth_diminfo @@ -614,7 +614,7 @@ subroutine specify_variable(this, name, dim_indices, global_level_data_size, loc character(len=*) , intent(in) :: name integer , intent(in) :: dim_indices(:) integer , intent(in) :: global_level_data_size - real(kind=8) , intent(inout), target :: local_data(:,:) ! todo: be able to set precision? + real(kind=WP) , intent(inout), target :: local_data(:,:) ! todo: be able to set precision? logical , intent(in) :: is_elem_based character(len=*) , intent(in) :: units, longname integer , intent(in) , optional :: ncat diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index 007cf2490..ea805d35a 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -3728,7 +3728,7 @@ subroutine io_r2g(n, partit, mesh) !___________________________________________________________________________ IF ((entry_x%accuracy == i_real8) .AND. (entry_y%accuracy == i_real8)) THEN -!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I, J, xmean, ymean) +!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I, J, temp_x, temp_y, xmean, ymean) DO J=1, size(entry_x%local_values_r8,dim=2) if (entry_x%is_elem_based) then xmean=sum(mesh%coord_nod2D(1, mesh%elem2D_nodes(:, J)))/3._WP @@ -3738,7 +3738,12 @@ subroutine io_r2g(n, partit, mesh) ymean=mesh%coord_nod2D(2, J) end if DO I=1, size(entry_x%local_values_r8,dim=1) - call vector_r2g(entry_x%local_values_r8(I,J), entry_y%local_values_r8(I,J), xmean, ymean, 0) + ! vector_r2g works in WP; round-trip through WP temps so the on-disk r8 buffer rotates correctly at WP=4 and WP=8 + temp_x=real(entry_x%local_values_r8(I,J), kind=WP) + temp_y=real(entry_y%local_values_r8(I,J), kind=WP) + call vector_r2g(temp_x, temp_y, xmean, ymean, 0) + entry_x%local_values_r8(I,J)=real(temp_x, real64) + entry_y%local_values_r8(I,J)=real(temp_y, real64) END DO END DO !$OMP END PARALLEL DO diff --git a/src/io_restart_file_group.F90 b/src/io_restart_file_group.F90 index 27e34a5f7..d6e340492 100644 --- a/src/io_restart_file_group.F90 +++ b/src/io_restart_file_group.F90 @@ -40,7 +40,7 @@ subroutine def_node_var_2d(this, name, longname, units, local_data, mesh, partit class(restart_file_group), target, intent(inout) :: this character(len=*), intent(in) :: name character(len=*), intent(in) :: units, longname - real(kind=8), target, intent(inout) :: local_data(:) ! todo: be able to set precision + real(kind=WP), target, intent(inout) :: local_data(:) ! todo: be able to set precision type(t_mesh), intent(in) :: mesh type(t_partit), intent(in) :: partit ! EO parameters @@ -54,7 +54,7 @@ subroutine def_node_var_3d(this, name, longname, units, local_data, mesh, partit class(restart_file_group), intent(inout) :: this character(len=*) , intent(in) :: name character(len=*) , intent(in) :: units, longname - real(kind=8) , intent(inout), target :: local_data(:,:) ! todo: be able to set precision + real(kind=WP) , intent(inout), target :: local_data(:,:) ! todo: be able to set precision type(t_mesh) , intent(in) :: mesh type(t_partit) , intent(in) :: partit integer , intent(in) , optional :: ncat @@ -78,7 +78,7 @@ subroutine def_elem_var_2d(this, name, longname, units, local_data, mesh, partit class(restart_file_group), intent(inout) :: this character(len=*), intent(in) :: name character(len=*), intent(in) :: units, longname - real(kind=8), target, intent(inout) :: local_data(:) ! todo: be able to set precision + real(kind=WP), target, intent(inout) :: local_data(:) ! todo: be able to set precision type(t_mesh), intent(in) :: mesh type(t_partit), intent(in) :: partit ! EO parameters @@ -93,7 +93,7 @@ subroutine def_elem_var_3d(this, name, longname, units, local_data, mesh, partit class(restart_file_group), intent(inout) :: this character(len=*), intent(in) :: name character(len=*), intent(in) :: units, longname - real(kind=8), target, intent(inout) :: local_data(:,:) ! todo: be able to set precision + real(kind=WP), target, intent(inout) :: local_data(:,:) ! todo: be able to set precision type(t_mesh), intent(in) :: mesh type(t_partit), intent(in) :: partit ! EO parameters @@ -143,7 +143,7 @@ subroutine def_node_var_2d_optional(this, name, longname, units, local_data, mes class(restart_file_group), target, intent(inout) :: this character(len=*), intent(in) :: name character(len=*), intent(in) :: units, longname - real(kind=8), target, intent(inout) :: local_data(:) ! todo: be able to set precision + real(kind=WP), target, intent(inout) :: local_data(:) ! todo: be able to set precision type(t_mesh), intent(in) :: mesh type(t_partit), intent(in) :: partit ! EO parameters @@ -157,7 +157,7 @@ subroutine def_node_var_3d_optional(this, name, longname, units, local_data, mes class(restart_file_group), intent(inout) :: this character(len=*) , intent(in) :: name character(len=*) , intent(in) :: units, longname - real(kind=8) , intent(inout), target :: local_data(:,:) ! todo: be able to set precision + real(kind=WP) , intent(inout), target :: local_data(:,:) ! todo: be able to set precision type(t_mesh) , intent(in) :: mesh type(t_partit) , intent(in) :: partit integer , intent(in) , optional :: ncat @@ -181,7 +181,7 @@ subroutine def_elem_var_2d_optional(this, name, longname, units, local_data, mes class(restart_file_group), intent(inout) :: this character(len=*), intent(in) :: name character(len=*), intent(in) :: units, longname - real(kind=8), target, intent(inout) :: local_data(:) ! todo: be able to set precision + real(kind=WP), target, intent(inout) :: local_data(:) ! todo: be able to set precision type(t_mesh), intent(in) :: mesh type(t_partit), intent(in) :: partit ! EO parameters @@ -196,7 +196,7 @@ subroutine def_elem_var_3d_optional(this, name, longname, units, local_data, mes class(restart_file_group), intent(inout) :: this character(len=*), intent(in) :: name character(len=*), intent(in) :: units, longname - real(kind=8), target, intent(inout) :: local_data(:,:) ! todo: be able to set precision + real(kind=WP), target, intent(inout) :: local_data(:,:) ! todo: be able to set precision type(t_mesh), intent(in) :: mesh type(t_partit), intent(in) :: partit ! EO parameters diff --git a/src/mod_transit.F90 b/src/mod_transit.F90 index d328fe4e8..329e29384 100644 --- a/src/mod_transit.F90 +++ b/src/mod_transit.F90 @@ -2,17 +2,18 @@ MODULE mod_transit ! Parameters, variables and functions for transient tracer simulations. ! By mbutzin, 2019-2021 + use o_PARAM, only: WP implicit none save ! Atmospheric pressure, local (dummy) variable, global-mean SLP, and local wind speed at 10 m - real(kind=8) :: press_a, mean_slp = 1.01325e5, wind_2 + real(kind=WP) :: press_a, mean_slp = 1.01325e5, wind_2 ! Atmospheric trace gas (dummy) values used in air-sea flux calculations ! Isotopic ratios are normalized and fractionation-corrected, ! volume mixing ratios are mole fractions in dry air. - real(kind=8) :: r14c_a = 1.0, & ! 14CO2 / 12CO2 (may vary with latitude in transient runs) + real(kind=WP) :: r14c_a = 1.0, & ! 14CO2 / 12CO2 (may vary with latitude in transient runs) r39ar_a = 1.0, & ! 39Ar / 40 Ar (homogeneous) xarg_a = 9.34e-3, & ! Argon (homogeneous) xCO2_a = 284.32e-6, & ! CO2 (CMIP6 & OMIP-BGC: 284.32e-6 for 1700-1850, PMIP4: 190.00e-6 for 21 ka BP) @@ -21,7 +22,7 @@ MODULE mod_transit xsf6_a = 0.0 ! SF6 (latitude dependent) ! Transient values of atmospheric trace gases (1d-arrays of variable length to be specified in namelist.config -> length_transit) - real(kind=8), allocatable, dimension(:) :: r14c_nh, r14c_tz, r14c_sh, & ! 14CO2 / 12CO2, latitude-dependent (e.g., bomb 14C) + real(kind=WP), allocatable, dimension(:) :: r14c_nh, r14c_tz, r14c_sh, & ! 14CO2 / 12CO2, latitude-dependent (e.g., bomb 14C) r14c_ti, & ! 14CO2 / 12CO2, homogenous (e.g., IntCal) xCO2_ti, & ! CO2 xf11_nh, xf11_sh, & ! CFC-11, latitude-dependent @@ -42,16 +43,16 @@ MODULE mod_transit ! Parameters which can be changed via namelist.oce (-> transit_param) ! Global-mean concentrations of DIC and Argon in the mixed layer (mol / m**3) - real(kind=8) :: dic_0 = 2.00, & ! GLODAPv2, 0-50 m: TCO2 ~ 2050 umol / kg + real(kind=WP) :: dic_0 = 2.00, & ! GLODAPv2, 0-50 m: TCO2 ~ 2050 umol / kg arg_0 = 0.01 ! Hamme et al. 2019, doi:10.1146/annurev-marine-121916-063604 ! Radioactive decay constants (1 / s; default values assume that 1 year = 365.00 days) - real(kind=8) :: decay14 = 3.8561e-12 , & ! 14C; t1/2 = 5700 a following OMIP-BGC + real(kind=WP) :: decay14 = 3.8561e-12 , & ! 14C; t1/2 = 5700 a following OMIP-BGC decay39 = 8.1708e-11 ! 39Ar; t1/2 = 269 a ! Further internal parameters ! Latitude of atmospheric boundary conditions and latitudinal interpolation weight - real(kind=8) :: y_abc, yy_nh + real(kind=WP) :: y_abc, yy_nh ! Tracer indices of transient tracers integer :: id_r14c, id_r39ar, id_f11, id_f12, id_sf6, index_transit_r14c, index_transit_r39ar, index_transit_f11, index_transit_f12, index_transit_sf6 ! Time index (=year) in transient simulations @@ -83,11 +84,11 @@ function iso_flux(which_gas, temp_c, sal, wind_2, f_ice, p_atm, x_gas, r_air, r_ ! for the abundant isotopologue. Positive values mean oceanic uptake. implicit none - real(kind=8) :: iso_flux + real(kind=WP) :: iso_flux ! Input parameters character(len=3), intent(in) :: which_gas ! trace gas name - real(kind=8), intent(in) :: temp_c, sal, & ! SST (deg C) and SSS ("PSU" or permil) + real(kind=WP), intent(in) :: temp_c, sal, & ! SST (deg C) and SSS ("PSU" or permil) wind_2, & ! wind speed at 10 m heigth squared f_ice, & ! sea-ice fractional coverage p_atm, & ! total atmospheric pressure (Pa) @@ -106,17 +107,17 @@ function gas_flux(which_gas, temp_c, sal, wind_2, f_ice, p_atm, x_gas, c_surf) ! Computes air-sea exchange gas fluxes in mol / (m**2 * s) , positive values mean oceanic uptake. implicit none - real(kind=8) :: gas_flux + real(kind=WP) :: gas_flux ! Input parameters character(len=3), intent(in) :: which_gas ! trace gas name - real(kind=8), intent(in) :: temp_c, sal, & ! SST (deg C) and SSS ("PSU" or permil) + real(kind=WP), intent(in) :: temp_c, sal, & ! SST (deg C) and SSS ("PSU" or permil) wind_2, & ! wind speed at 10 m heigth squared f_ice, & ! sea-ice fractional coverage p_atm, & ! total atmospheric pressure (Pa) x_gas, & ! atmospheric mole fraction c_surf ! marine surface water concentration (mol / m**3) ! Internal variables - real(kind=8) :: c_sat ! marine saturation concentration (mol / m**3) + real(kind=WP) :: c_sat ! marine saturation concentration (mol / m**3) c_sat = solub(which_gas, temp_c, sal) * p_atm / 1.01325e5 * x_gas gas_flux = transfer_vel(which_gas, temp_c, wind_2) * (c_sat - c_surf) * (1. - f_ice) @@ -128,12 +129,12 @@ function solub(which_gas, temp_c, sal) ! Computes the solubility of trace gases in seawater. ! This parametrization includes the effect of water vapor. implicit none - real(kind=8) :: solub ! solubility ((p)mol / (m**3 * atm)) + real(kind=WP) :: solub ! solubility ((p)mol / (m**3 * atm)) ! Input parameters character(len=3), intent(in) :: which_gas ! tracer name - real(kind=8), intent(in) :: temp_c, & ! temperature (deg C) + real(kind=WP), intent(in) :: temp_c, & ! temperature (deg C) sal ! salinity ("PSU" or permil) - real(kind=8) :: a1, a2, a3, a4, & ! polynomial coefficients of the + real(kind=WP) :: a1, a2, a3, a4, & ! polynomial coefficients of the b1, b2, b3, b4, c1, & ! solubility function temp_k100, & ! water temperature in K / 100 con2con ! concentration units conversion factor @@ -182,12 +183,12 @@ function sc_660(which_gas, temp_c) ! normalized to 20 degC (Sc(CO2) ~660; Wanninkhof 2014, tab. 1)). implicit none ! Result - real(kind=8) :: sc_660 ! Schmidt number + real(kind=WP) :: sc_660 ! Schmidt number ! Input parameters character(len=3), intent(in) :: which_gas ! tracer name - real(kind=8), intent(in) :: temp_c ! temperature (deg C) + real(kind=WP), intent(in) :: temp_c ! temperature (deg C) ! Internal parameters and/or variables - real(kind=8) :: as, bs, cs, ds, es ! polynomial coefficients + real(kind=WP) :: as, bs, cs, ds, es ! polynomial coefficients select case (which_gas) case ("co2") ! CO2 @@ -212,10 +213,10 @@ function transfer_vel(which_gas, temp_c, wind_2) ! Compute gas transfer velocities of / for tracers implicit none ! Result - real(kind=8) :: transfer_vel ! transfer velocity (m / s) + real(kind=WP) :: transfer_vel ! transfer velocity (m / s) ! Input parameters character(len=3), intent(in) :: which_gas ! tracer name - real(kind=8), intent(in) :: temp_c, & ! temperature (deg C) + real(kind=WP), intent(in) :: temp_c, & ! temperature (deg C) wind_2 ! wind speed squared at 10 m height (m / s) ! Wanninkhof (2014), eq. (4) with a = 0.251 (cm / h) / (m / s)**2 -> 6.9722e-7 s / m @@ -232,20 +233,20 @@ function speed_2(windstr_x, windstr_y) ! We follow Peixoto & Oort (1992, Eq. (10.28), (10,29)) and Charnock (1955); ! also see MPI report 349 (2003), Eq. (5.7). implicit none - real(kind=8) :: speed_2 + real(kind=WP) :: speed_2 ! Input - real(kind=8), intent(in) :: windstr_x, windstr_y + real(kind=WP), intent(in) :: windstr_x, windstr_y ! Internal variables and parameters ! Zonal and meridional velocities at 10 m height - real(kind=8) :: u_10, v_10 + real(kind=WP) :: u_10, v_10 ! Zonal and meridional friction velocities - real(kind=8) :: u_fric, v_fric + real(kind=WP) :: u_fric, v_fric ! Zonal and meridional roughness lengths - real(kind=8) :: l_rough_x, l_rough_y + real(kind=WP) :: l_rough_x, l_rough_y ! Inverse von-Karman constant (0.4), Charnock constant (0.018) divided by g, inverse density of air (1.3), log(10) - real(kind=8), parameter :: inv_karm = 2.5, charn_g = 0.00173, inv_dens_air = 0.76923, log_10 = 2.30258 + real(kind=WP), parameter :: inv_karm = 2.5, charn_g = 0.00173, inv_dens_air = 0.76923, log_10 = 2.30258 ! Calculate friction velocities (Peixoto & Oort, 1992, Eq. (10.28)) u_fric = sqrt(abs(windstr_x) * inv_dens_air) @@ -272,7 +273,7 @@ subroutine read_transit_input(ifileunit) ! Internal variables integer, intent(in) :: ifileunit integer :: jj - real(kind=8), allocatable, dimension(:) :: d14c_nh, d14c_tz, d14c_sh, d14c_ti, d13c_dummy + real(kind=WP), allocatable, dimension(:) :: d14c_nh, d14c_tz, d14c_sh, d14c_ti, d13c_dummy if (anthro_transit) then ! Anthropogenic input for 1850 - 2023 CE diff --git a/src/oce_ale_pressure_bv.F90 b/src/oce_ale_pressure_bv.F90 index a6192ec20..2f02e4ebe 100644 --- a/src/oce_ale_pressure_bv.F90 +++ b/src/oce_ale_pressure_bv.F90 @@ -3323,9 +3323,9 @@ subroutine init_ref_density_advanced(tracers, partit, mesh) type(t_tracer), intent(in), target :: tracers integer :: node, nz, nzmin, nzmax real(kind=WP) :: rhopot, bulk_0, bulk_pz, bulk_pz2, rho - real(kind=8) :: T, S, auxz, x, y + real(kind=WP) :: T, S, auxz, x, y real(kind=WP), dimension(:,:), pointer :: temp, salt - real(kind=8) :: ref_temp1D(mesh%nl-1), ref_salt1D(mesh%nl-1), vol1D(mesh%nl-1) + real(kind=WP) :: ref_temp1D(mesh%nl-1), ref_salt1D(mesh%nl-1), vol1D(mesh%nl-1) #include "associate_part_def.h" #include "associate_mesh_def.h" #include "associate_part_ass.h" @@ -3351,9 +3351,9 @@ subroutine init_ref_density_advanced(tracers, partit, mesh) end do end do !$OMP END PARALLEL DO -call MPI_Allreduce(MPI_IN_PLACE, ref_temp1D, mesh%nl-1, MPI_DOUBLE, MPI_SUM, partit%MPI_COMM_FESOM, MPIerr) -call MPI_Allreduce(MPI_IN_PLACE, ref_salt1D, mesh%nl-1, MPI_DOUBLE, MPI_SUM, partit%MPI_COMM_FESOM, MPIerr) -call MPI_Allreduce(MPI_IN_PLACE, vol1D, mesh%nl-1, MPI_DOUBLE, MPI_SUM, partit%MPI_COMM_FESOM, MPIerr) +call MPI_Allreduce(MPI_IN_PLACE, ref_temp1D, mesh%nl-1, MPI_WP, MPI_SUM, partit%MPI_COMM_FESOM, MPIerr) +call MPI_Allreduce(MPI_IN_PLACE, ref_salt1D, mesh%nl-1, MPI_WP, MPI_SUM, partit%MPI_COMM_FESOM, MPIerr) +call MPI_Allreduce(MPI_IN_PLACE, vol1D, mesh%nl-1, MPI_WP, MPI_SUM, partit%MPI_COMM_FESOM, MPIerr) where( vol1D > 1.e-12_WP) !more than 0.! ref_temp1D=ref_temp1D/vol1D diff --git a/src/oce_dyn.F90 b/src/oce_dyn.F90 index d876bc95f..fa35c8a46 100755 --- a/src/oce_dyn.F90 +++ b/src/oce_dyn.F90 @@ -1003,8 +1003,8 @@ subroutine check_validviscopt_5(partit, mesh) !___________________________________________________________________________ ! compute global mean ratio --> core2 Ratio=4.26 (eddy parameterizted), ! dart Ratio=0.97 (eddy resolving/permitting) - call MPI_AllREDUCE(loc_R, glb_R, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) - call MPI_AllREDUCE(loc_A, glb_A, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc_R, glb_R, 1, MPI_WP, MPI_SUM, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc_A, glb_A, 1, MPI_WP, MPI_SUM, MPI_COMM_FESOM, MPIerr) glb_R = glb_R/glb_A !___________________________________________________________________________ diff --git a/src/oce_mesh.F90 b/src/oce_mesh.F90 index 03f2f41d8..3ab39f419 100755 --- a/src/oce_mesh.F90 +++ b/src/oce_mesh.F90 @@ -403,8 +403,8 @@ SUBROUTINE read_mesh(partit, mesh) end if end do end if - call MPI_BCast(rbuff(1:k,1), k, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) - call MPI_BCast(rbuff(1:k,2), k, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(rbuff(1:k,1), k, MPI_WP, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(rbuff(1:k,2), k, MPI_WP, 0, MPI_COMM_FESOM, ierror) call MPI_BCast(ibuff(1:k,2), k, MPI_INTEGER, 0, MPI_COMM_FESOM, ierror) ! fill the local arrays do n=1, k @@ -574,7 +574,7 @@ SUBROUTINE read_mesh(partit, mesh) end if allocate(mesh%zbar(mesh%nl)) ! allocate the array for storing the standard depths if (mype==0) read(fileID,*) mesh%zbar - call MPI_BCast(mesh%zbar, mesh%nl, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(mesh%zbar, mesh%nl, MPI_WP, 0, MPI_COMM_FESOM, ierror) if(mesh%zbar(2)>0) mesh%zbar=-mesh%zbar ! zbar is negative allocate(mesh%Z(mesh%nl-1)) mesh%Z=mesh%zbar(1:mesh%nl-1)+mesh%zbar(2:mesh%nl) ! mid-depths of cells @@ -615,7 +615,7 @@ SUBROUTINE read_mesh(partit, mesh) end if allocate(mesh%zbar(mesh%nl)) ! allocate the array for storing the standard depths if (mype==0) read(fileID,*) mesh%zbar - call MPI_BCast(mesh%zbar, mesh%nl, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(mesh%zbar, mesh%nl, MPI_WP, 0, MPI_COMM_FESOM, ierror) if(mesh%zbar(2)>0) mesh%zbar=-mesh%zbar ! zbar is negative allocate(mesh%Z(mesh%nl-1)) mesh%Z=mesh%zbar(1:mesh%nl-1)+mesh%zbar(2:mesh%nl) ! mid-depths of cells @@ -693,7 +693,7 @@ SUBROUTINE read_mesh(partit, mesh) ! the maximum depth on earth in marianen trench if ( flag_wrongaux3d==0 .and. any(abs(rbuff(1:k,1))>11000.0_WP) ) flag_wrongaux3d=1 end if - call MPI_BCast(rbuff(1:k,1), k, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(rbuff(1:k,1), k, MPI_WP, 0, MPI_COMM_FESOM, ierror) do n=1, k x=rbuff(n,1) @@ -746,7 +746,7 @@ SUBROUTINE read_mesh(partit, mesh) ! the maximum depth on earth in marianen trench if ( flag_wrongaux3d==0 .and. any(abs(rbuff(1:k,1))>11000.0_WP) ) flag_wrongaux3d=1 end if - call MPI_BCast(rbuff(1:k,1), k, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(rbuff(1:k,1), k, MPI_WP, 0, MPI_COMM_FESOM, ierror) do n=1, k x=rbuff(n,1) @@ -1470,7 +1470,7 @@ subroutine find_levels_cavity(partit, mesh) !___________________________________________________________________ ! broadcast chunk buffer to all other CPUs (k...size of buffer) - call MPI_BCast(rbuff(1:k), k, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(rbuff(1:k), k, MPI_WP, 0, MPI_COMM_FESOM, ierror) !___________________________________________________________________ ! fill the local arrays @@ -1543,7 +1543,7 @@ subroutine find_levels_cavity(partit, mesh) !___________________________________________________________________ ! broadcast chunk buffer to all other CPUs (k...size of buffer) - call MPI_BCast(rbuff(1:k), k, MPI_DOUBLE_PRECISION, 0, MPI_COMM_FESOM, ierror) + call MPI_BCast(rbuff(1:k), k, MPI_WP, 0, MPI_COMM_FESOM, ierror) !___________________________________________________________________ ! fill the local arrays @@ -2395,9 +2395,9 @@ SUBROUTINE mesh_areas(partit, mesh) end do mesh%ocean_area=0.0_WP mesh%ocean_areawithcav=0.0_WP - call MPI_AllREDUCE(vol, mesh%ocean_area, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & + call MPI_AllREDUCE(vol, mesh%ocean_area, 1, MPI_WP, MPI_SUM, & MPI_COMM_FESOM, MPIerr) - call MPI_AllREDUCE(vol2, mesh%ocean_areawithcav, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & + call MPI_AllREDUCE(vol2, mesh%ocean_areawithcav, 1, MPI_WP, MPI_SUM, & MPI_COMM_FESOM, MPIerr) !___write mesh statistics___________________________________________________ @@ -2821,7 +2821,7 @@ SUBROUTINE check_mesh_consistency(partit, mesh) aux(nz)=aux(nz)+mesh%areasvol(nz, n) end do end do - call MPI_AllREDUCE(aux, vol_n, mesh%nl, MPI_DOUBLE_PRECISION, MPI_SUM, & + call MPI_AllREDUCE(aux, vol_n, mesh%nl, MPI_WP, MPI_SUM, & MPI_COMM_FESOM, MPIerr) aux=0._WP @@ -2832,7 +2832,7 @@ SUBROUTINE check_mesh_consistency(partit, mesh) aux(nz)=aux(nz)+mesh%elem_area(elem) end do end do - call MPI_AllREDUCE(aux, vol_e, mesh%nl, MPI_DOUBLE_PRECISION, MPI_SUM, & + call MPI_AllREDUCE(aux, vol_e, mesh%nl, MPI_WP, MPI_SUM, & MPI_COMM_FESOM, MPIerr) if (mype==0) then @@ -2879,7 +2879,7 @@ subroutine check_total_volume(partit, mesh) aux=aux+areasvol(nz, n)*hnode(nz,n) end do end do - call MPI_AllREDUCE(aux, vol_n, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(aux, vol_n, 1, MPI_WP, MPI_SUM, MPI_COMM_FESOM, MPIerr) !___________________________________________________________________________ ! total ocean volume on elements aux=0._WP @@ -2890,7 +2890,7 @@ subroutine check_total_volume(partit, mesh) aux=aux+elem_area(elem)*helem(nz,elem) end do end do - call MPI_AllREDUCE(aux, vol_e, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(aux, vol_e, 1, MPI_WP, MPI_SUM, MPI_COMM_FESOM, MPIerr) !___write mesh statistics___________________________________________________ if (mype==0) then diff --git a/src/oce_modules.F90 b/src/oce_modules.F90 index c80a1d1e3..dfa0f8cfe 100755 --- a/src/oce_modules.F90 +++ b/src/oce_modules.F90 @@ -5,7 +5,12 @@ !========================================================== MODULE o_PARAM -integer, parameter :: WP=8 ! Working precision +! Working precision - configurable via preprocessor (-DUSE_SINGLE_PRECISION) +#if defined(USE_SINGLE_PRECISION) +integer, parameter :: WP=4 ! Single precision +#else +integer, parameter :: WP=8 ! Double precision (default) +#endif integer, parameter :: MAX_PATH=4096 ! Maximum file path length integer :: mstep real(kind=WP), parameter :: pi=3.14159265358979 diff --git a/src/solver.F90 b/src/solver.F90 index c70876e23..dbd4f1176 100644 --- a/src/solver.F90 +++ b/src/solver.F90 @@ -150,7 +150,7 @@ subroutine ssh_solve_cg(x, rhs, solverinfo, partit, mesh) s_old = sum(rhs(1:myDim_nod2D) * rhs(1:myDim_nod2D)) #endif - call MPI_Allreduce(MPI_IN_PLACE, s_old, 1, MPI_DOUBLE, MPI_SUM, partit%MPI_COMM_FESOM, MPIerr) + call MPI_Allreduce(MPI_IN_PLACE, s_old, 1, MPI_WP, MPI_SUM, partit%MPI_COMM_FESOM, MPIerr) rtol=solverinfo%soltol*sqrt(s_old/real(nod2D,WP)) ! ============== ! Compute r0 @@ -188,7 +188,7 @@ subroutine ssh_solve_cg(x, rhs, solverinfo, partit, mesh) s_old = sum(rr(1:myDim_nod2D) * zz(1:myDim_nod2D)) #endif - call MPI_Allreduce(MPI_IN_PLACE, s_old, 1, MPI_DOUBLE, MPI_SUM, partit%MPI_COMM_FESOM, MPIerr) + call MPI_Allreduce(MPI_IN_PLACE, s_old, 1, MPI_WP, MPI_SUM, partit%MPI_COMM_FESOM, MPIerr) ! =============== ! Iterations @@ -218,7 +218,7 @@ subroutine ssh_solve_cg(x, rhs, solverinfo, partit, mesh) s_aux = sum(pp(1:myDim_nod2D) * App(1:myDim_nod2D)) #endif - call MPI_Allreduce(MPI_IN_PLACE, s_aux, 1, MPI_DOUBLE, MPI_SUM, partit%MPI_COMM_FESOM, MPIerr) + call MPI_Allreduce(MPI_IN_PLACE, s_aux, 1, MPI_WP, MPI_SUM, partit%MPI_COMM_FESOM, MPIerr) al=s_old/s_aux ! =========== @@ -256,7 +256,7 @@ subroutine ssh_solve_cg(x, rhs, solverinfo, partit, mesh) sprod(1) = sum(rr(1:myDim_nod2D) * rr(1:myDim_nod2D)) #endif - call MPI_Allreduce(MPI_IN_PLACE, sprod, 2, MPI_DOUBLE, MPI_SUM, partit%MPI_COMM_FESOM, MPIerr) + call MPI_Allreduce(MPI_IN_PLACE, sprod, 2, MPI_WP, MPI_SUM, partit%MPI_COMM_FESOM, MPIerr) !$OMP BARRIER ! =========== diff --git a/src/write_step_info.F90 b/src/write_step_info.F90 index 0ad26ee50..d05362f3b 100644 --- a/src/write_step_info.F90 +++ b/src/write_step_info.F90 @@ -133,11 +133,11 @@ subroutine write_step_info(istep, outfreq, ice, dynamics, tracers, partit, mesh) end if !_______________________________________________________________________ - call MPI_AllREDUCE(loc_eta , int_eta , 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) - call MPI_AllREDUCE(loc_hbar , int_hbar , 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) -!PS call MPI_AllREDUCE(loc_deta , int_deta , 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) - call MPI_AllREDUCE(loc_dhbar, int_dhbar, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) - call MPI_AllREDUCE(loc_wflux, int_wflux, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc_eta , int_eta , 1, MPI_WP, MPI_SUM, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc_hbar , int_hbar , 1, MPI_WP, MPI_SUM, MPI_COMM_FESOM, MPIerr) +!PS call MPI_AllREDUCE(loc_deta , int_deta , 1, MPI_WP, MPI_SUM, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc_dhbar, int_dhbar, 1, MPI_WP, MPI_SUM, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc_wflux, int_wflux, 1, MPI_WP, MPI_SUM, MPI_COMM_FESOM, MPIerr) int_eta = int_eta /ocean_areawithcav int_hbar = int_hbar /ocean_areawithcav @@ -146,89 +146,89 @@ subroutine write_step_info(istep, outfreq, ice, dynamics, tracers, partit, mesh) int_wflux= int_wflux/ocean_areawithcav !_______________________________________________________________________ loc=omp_min_max_sum1(eta_n, 1, myDim_nod2D, 'min', partit) - call MPI_AllREDUCE(loc , min_eta , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , min_eta , 1, MPI_WP, MPI_MIN, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(hbar, 1, myDim_nod2D, 'min', partit) - call MPI_AllREDUCE(loc , min_hbar , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , min_hbar , 1, MPI_WP, MPI_MIN, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(water_flux, 1, myDim_nod2D, 'min', partit) - call MPI_AllREDUCE(loc , min_wflux, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , min_wflux, 1, MPI_WP, MPI_MIN, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(heat_flux, 1, myDim_nod2D, 'min', partit) - call MPI_AllREDUCE(loc , min_hflux, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , min_hflux, 1, MPI_WP, MPI_MIN, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum2(tracers%data(1)%values, 1, nl-1, 1, myDim_nod2D, 'min', partit, 0.0_WP) - call MPI_AllREDUCE(loc , min_temp , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , min_temp , 1, MPI_WP, MPI_MIN, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum2(tracers%data(2)%values, 1, nl-1, 1, myDim_nod2D, 'min', partit, 0.0_WP) - call MPI_AllREDUCE(loc , min_salt , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , min_salt , 1, MPI_WP, MPI_MIN, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(Wvel(1,:), 1, myDim_nod2D, 'min', partit) - call MPI_AllREDUCE(loc , min_wvel , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , min_wvel , 1, MPI_WP, MPI_MIN, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(Wvel(2,:), 1, myDim_nod2D, 'min', partit) - call MPI_AllREDUCE(loc , min_wvel2 , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , min_wvel2 , 1, MPI_WP, MPI_MIN, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(UVnode(1,1,:), 1, myDim_nod2D, 'min', partit) - call MPI_AllREDUCE(loc , min_uvel , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , min_uvel , 1, MPI_WP, MPI_MIN, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(UVnode(1,2,:), 1, myDim_nod2D, 'min', partit) - call MPI_AllREDUCE(loc , min_uvel2, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , min_uvel2, 1, MPI_WP, MPI_MIN, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(UVnode(2,1,:), 1, myDim_nod2D, 'min', partit) - call MPI_AllREDUCE(loc , min_vvel , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , min_vvel , 1, MPI_WP, MPI_MIN, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(UVnode(2,2,:), 1, myDim_nod2D, 'min', partit) - call MPI_AllREDUCE(loc , min_vvel2 , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , min_vvel2 , 1, MPI_WP, MPI_MIN, MPI_COMM_FESOM, MPIerr) if ( .not. dynamics%use_ssh_se_subcycl) then loc=omp_min_max_sum1(d_eta, 1, myDim_nod2D, 'min', partit) else loc=omp_min_max_sum1(hbar-hbar_old, 1, myDim_nod2D, 'min', partit) end if - call MPI_AllREDUCE(loc , min_deta , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , min_deta , 1, MPI_WP, MPI_MIN, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(hnode(1,:), 1, myDim_nod2D, 'min', partit) - call MPI_AllREDUCE(loc , min_hnode , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , min_hnode , 1, MPI_WP, MPI_MIN, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(hnode(2,:), 1, myDim_nod2D, 'min', partit) - call MPI_AllREDUCE(loc , min_hnode2 , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , min_hnode2 , 1, MPI_WP, MPI_MIN, MPI_COMM_FESOM, MPIerr) !_______________________________________________________________________ loc=omp_min_max_sum1(eta_n, 1, myDim_nod2D, 'max', partit) - call MPI_AllREDUCE(loc , max_eta , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_eta , 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(hbar, 1, myDim_nod2D, 'max', partit) - call MPI_AllREDUCE(loc , max_hbar , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_hbar , 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(water_flux, 1, myDim_nod2D, 'max', partit) - call MPI_AllREDUCE(loc , max_wflux, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_wflux, 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(heat_flux, 1, myDim_nod2D, 'max', partit) - call MPI_AllREDUCE(loc , max_hflux, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_hflux, 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum2(tracers%data(1)%values, 1, nl-1, 1, myDim_nod2D, 'max', partit, 0.0_WP) - call MPI_AllREDUCE(loc , max_temp , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_temp , 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum2(tracers%data(2)%values, 1, nl-1, 1, myDim_nod2D, 'max', partit, 0.0_WP) - call MPI_AllREDUCE(loc , max_salt , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_salt , 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(Wvel(1,:), 1, myDim_nod2D, 'max', partit) - call MPI_AllREDUCE(loc , max_wvel , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_wvel , 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(Wvel(2,:), 1, myDim_nod2D, 'max', partit) - call MPI_AllREDUCE(loc , max_wvel2 , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_wvel2 , 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(UVnode(1,1,:), 1, myDim_nod2D, 'max', partit) - call MPI_AllREDUCE(loc , max_uvel , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_uvel , 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(UVnode(1,2,:), 1, myDim_nod2D, 'max', partit) - call MPI_AllREDUCE(loc , max_uvel2 , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_uvel2 , 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(UVnode(2,1,:), 1, myDim_nod2D, 'max', partit) - call MPI_AllREDUCE(loc , max_vvel , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_vvel , 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(UVnode(2,2,:), 1, myDim_nod2D, 'max', partit) - call MPI_AllREDUCE(loc , max_vvel2 , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_vvel2 , 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) if ( .not. dynamics%use_ssh_se_subcycl) then loc=omp_min_max_sum1(d_eta, 1, myDim_nod2D, 'max', partit) else loc=omp_min_max_sum1(hbar-hbar_old, 1, myDim_nod2D, 'max', partit) end if - call MPI_AllREDUCE(loc , max_deta , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_deta , 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(hnode(1, :), 1, myDim_nod2D, 'max', partit) - call MPI_AllREDUCE(loc , max_hnode , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_hnode , 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum1(hnode(2, :), 1, myDim_nod2D, 'max', partit) - call MPI_AllREDUCE(loc , max_hnode2 , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_hnode2 , 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum2(CFL_z, 1, nl-1, 1, myDim_nod2D, 'max', partit) - call MPI_AllREDUCE(loc , max_cfl_z , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_cfl_z , 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum2(pgf_x, 1, nl-1, 1, myDim_nod2D, 'max', partit) - call MPI_AllREDUCE(loc , max_pgfx , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_pgfx , 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum2(pgf_y, 1, nl-1, 1, myDim_nod2D, 'max', partit) - call MPI_AllREDUCE(loc , max_pgfy , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_pgfy , 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) if (use_ice) then loc=omp_min_max_sum1(m_ice, 1, myDim_nod2D, 'max', partit) - call MPI_AllREDUCE(loc , max_m_ice , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_m_ice , 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) end if loc=omp_min_max_sum2(Av, 1, nl, 1, myDim_elem2D, 'max', partit) - call MPI_AllREDUCE(loc , max_av , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_av , 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) loc=omp_min_max_sum2(Kv, 1, nl, 1, myDim_nod2D, 'max', partit) - call MPI_AllREDUCE(loc , max_kv , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + call MPI_AllREDUCE(loc , max_kv , 1, MPI_WP, MPI_MAX, MPI_COMM_FESOM, MPIerr) !_______________________________________________________________________ if (mype==0) then write(*,*) '___CHECK GLOBAL OCEAN VARIABLES --> mstep=',mstep From 48c373285c8c7d8e111cf799cedf29075fdfdc9b Mon Sep 17 00:00:00 2001 From: Suvarchal Kumar Cheedela Date: Tue, 30 Jun 2026 05:06:00 +0200 Subject: [PATCH 2/8] oce: fix single-precision NaN in KPP vertical mixing (epsln underflow) The KPP divide-by-zero guard epsln=1.0e-40 is a denormal in real32 and is flushed to zero under Intel -O3 (FTZ/DAZ) in single-precision builds, defeating the /(wm+epsln), /(hbl+epsln), /(dVsq+Vtsq+epsln) guards and producing 0/0 NaN at the cold start (Wvel NaN at step 1 on core2). Make epsln precision-dependent: 1.0e-20 (normal, FTZ-safe) in single precision, 1.0e-40 (unchanged) in double precision. The double-precision build is bit-identical. --- src/oce_ale_mixing_kpp.F90 | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/oce_ale_mixing_kpp.F90 b/src/oce_ale_mixing_kpp.F90 index c7204cd46..2da1bb995 100755 --- a/src/oce_ale_mixing_kpp.F90 +++ b/src/oce_ale_mixing_kpp.F90 @@ -47,7 +47,17 @@ MODULE o_mixing_KPP_mod real(KIND=WP), dimension(:,:), allocatable :: ghats ! nonlocal transport (s/m^2) real(KIND=WP), dimension(:), allocatable :: hbl ! boundary layer depth + ! KPP divide-by-zero guard added to denominators (/(wm+epsln), /(hbl+epsln), + ! /(dVsq+Vtsq+epsln), (ws*hbl+epsln), ...). It MUST stay a NORMAL number in the + ! active working precision: 1.0e-40 is a denormal in real32 and Intel -O3 flushes + ! it to 0 (FTZ/DAZ), which defeats every guard -> 0/0 = NaN in single precision at + ! the cold start. Single precision therefore uses a value that survives FTZ; double + ! precision keeps 1.0e-40 so the default build is numerically unchanged. +#if defined(USE_SINGLE_PRECISION) + real(KIND=WP), parameter :: epsln = 1.0e-20_WP ! a small value (SP: normal, FTZ-safe) +#else real(KIND=WP), parameter :: epsln = 1.0e-40_WP ! a small value +#endif real(KIND=WP), parameter :: epsilon_kpp = 0.1_WP real(KIND=WP), parameter :: vonk = 0.4_WP From e88de1acb42f390f3fc42f3412cb9c8348a45b58 Mon Sep 17 00:00:00 2001 From: Suvarchal Kumar Cheedela Date: Tue, 30 Jun 2026 05:06:43 +0200 Subject: [PATCH 3/8] testing: raise stack limit when launching fesom.x / fesom_meshdiag Intel-compiled fesom.x stack-allocates large per-column automatic arrays and SIGSEGVs on global meshes (e.g. core2) under the default 8 MB stack limit. Wrap the test launch in /bin/sh -c "ulimit -s unlimited || ulimit -s 1048576; exec ..." so mpiexec and the ranks it spawns inherit a larger stack. Degrades gracefully where unlimited is disallowed; harmless for GNU/CI. --- cmake/FesomTesting.cmake | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/cmake/FesomTesting.cmake b/cmake/FesomTesting.cmake index f99f97b9f..e2e72fb83 100644 --- a/cmake/FesomTesting.cmake +++ b/cmake/FesomTesting.cmake @@ -409,9 +409,13 @@ function(add_fesom_test_with_options TEST_NAME MESH_NAME STEP_PER_DAY RUN_LENGTH set(ENV{OMPI_MCA_rmaps_base_oversubscribe} \"1\") set(ENV{PRTE_MCA_rmaps_default_mapping_policy} \":oversubscribe\") - # Run FESOM with MPI + # Run FESOM with MPI. + # Wrap the launch in a shell that raises the stack limit: Intel-compiled + # fesom.x puts large per-column automatic arrays on the stack and SIGSEGVs + # on global meshes (e.g. core2) under the default 8 MB limit. ulimit is + # inherited by mpiexec and the ranks it spawns locally. Harmless for GNU/CI. execute_process( - COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} ${FESOM_TEST_NP} ${CMAKE_BINARY_DIR}/bin/fesom.x + COMMAND /bin/sh -c \"ulimit -s unlimited 2>/dev/null || ulimit -s 1048576 2>/dev/null; exec '${MPIEXEC_EXECUTABLE}' ${MPIEXEC_NUMPROC_FLAG} ${FESOM_TEST_NP} '${CMAKE_BINARY_DIR}/bin/fesom.x'\" WORKING_DIRECTORY \"${TEST_RUN_DIR}\" RESULT_VARIABLE test_result OUTPUT_VARIABLE test_output @@ -441,9 +445,10 @@ function(add_fesom_test_with_options TEST_NAME MESH_NAME STEP_PER_DAY RUN_LENGTH file(MAKE_DIRECTORY \"${TEST_RUN_DIR}\") file(MAKE_DIRECTORY \"${RESULT_DIR}\") - # Run FESOM + # Run FESOM (serial). Raise the stack limit as for the MPI case above + # (Intel fesom.x overflows the default 8 MB stack on large meshes). execute_process( - COMMAND ${CMAKE_BINARY_DIR}/bin/fesom.x + COMMAND /bin/sh -c \"ulimit -s unlimited 2>/dev/null || ulimit -s 1048576 2>/dev/null; exec '${CMAKE_BINARY_DIR}/bin/fesom.x'\" WORKING_DIRECTORY \"${TEST_RUN_DIR}\" RESULT_VARIABLE test_result OUTPUT_VARIABLE test_output @@ -584,9 +589,10 @@ function(add_fesom_meshdiag_test_with_options TEST_NAME MESH_NAME RUNID) set(ENV{OMPI_MCA_rmaps_base_oversubscribe} \"1\") set(ENV{PRTE_MCA_rmaps_default_mapping_policy} \":oversubscribe\") - # Run fesom_meshdiag with MPI + # Run fesom_meshdiag with MPI. Raise the stack limit (Intel builds overflow + # the default 8 MB stack on large meshes); inherited by mpiexec + ranks. execute_process( - COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} ${FESOM_TEST_NP} ${CMAKE_BINARY_DIR}/bin/fesom_meshdiag + COMMAND /bin/sh -c \"ulimit -s unlimited 2>/dev/null || ulimit -s 1048576 2>/dev/null; exec '${MPIEXEC_EXECUTABLE}' ${MPIEXEC_NUMPROC_FLAG} ${FESOM_TEST_NP} '${CMAKE_BINARY_DIR}/bin/fesom_meshdiag'\" WORKING_DIRECTORY \"${TEST_RUN_DIR}\" RESULT_VARIABLE test_result OUTPUT_VARIABLE test_output From 01fa1cfd3990ae2506a371d8b8d5f0985cd57ede Mon Sep 17 00:00:00 2001 From: Suvarchal Kumar Cheedela Date: Tue, 30 Jun 2026 05:08:39 +0200 Subject: [PATCH 4/8] cvmix: run CVMix in double precision to support single-precision builds CVMix is a fixed double-precision library (kind cvmix_r8 = real64); it was forced OFF under USE_SINGLE_PRECISION. Link it in both precisions and convert WP<->cvmix_r8 at the CVMix call boundary in the wrappers (scalar args cast, array args via cvmix_r8 temp buffers; FESOM-side arrays stay WP). CVMix always runs in double precision, so the mixing-scheme menu is identical in single and double precision (enabling like-for-like comparison) and the default double-precision build is numerically unchanged. Verified on core2: KPP, PP, cvmix_PP, cvmix_KPP and cvmix_TKE all run in GNU and Intel single precision. --- CMakeLists.txt | 12 +- src/cvmix_driver/gen_modules_cvmix_idemix.F90 | 103 +++++-- src/cvmix_driver/gen_modules_cvmix_kpp.F90 | 266 ++++++++++++------ src/cvmix_driver/gen_modules_cvmix_pp.F90 | 42 +-- src/cvmix_driver/gen_modules_cvmix_tidal.F90 | 56 ++-- src/cvmix_driver/gen_modules_cvmix_tke.F90 | 139 ++++++--- 6 files changed, 416 insertions(+), 202 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index b7b5b6414..a5135434f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -27,12 +27,12 @@ set(VERBOSE OFF CACHE BOOL "toggle debug output") set(FESOM_PROFILING OFF CACHE BOOL "Enable enhanced profiling system") set(CVMIX ON CACHE BOOL "Download/Compile/Install/Link CVMix libray") -# Single precision (WP=4) build: CVMix is a double-precision library and is out of scope -# for the SP port, so disable it to keep the SP build clean (native KPP/PP still work). -if(USE_SINGLE_PRECISION) - message(STATUS "USE_SINGLE_PRECISION=ON: forcing CVMIX=OFF (out of scope for single precision)") - set(CVMIX OFF CACHE BOOL "Download/Compile/Install/Link CVMix libray" FORCE) -endif() +# CVMix is a fixed double-precision library (its kind is cvmix_r8 = real64). It is +# linked in BOTH double- and single-precision FESOM builds: the FESOM wrappers in +# src/cvmix_driver/ convert WP<->cvmix_r8 at the CVMix call boundary, so CVMix always +# runs in double precision. This keeps the available mixing schemes identical in SP and +# DP (so SP and DP configurations can be compared like-for-like) and leaves the default +# DP build numerically unchanged. CVMIX can still be turned off explicitly (-DCVMIX=OFF). # Testing options option(BUILD_TESTING "Build tests" OFF) diff --git a/src/cvmix_driver/gen_modules_cvmix_idemix.F90 b/src/cvmix_driver/gen_modules_cvmix_idemix.F90 index 615d1ae51..5125b5258 100644 --- a/src/cvmix_driver/gen_modules_cvmix_idemix.F90 +++ b/src/cvmix_driver/gen_modules_cvmix_idemix.F90 @@ -18,9 +18,14 @@ module g_cvmix_idemix !___________________________________________________________________________ ! module calls from cvmix library - use cvmix_idemix, only : init_idemix, calc_idemix_v0, cvmix_coeffs_idemix + use cvmix_idemix, only : init_idemix, calc_idemix_v0, cvmix_coeffs_idemix use cvmix_put_get, only : cvmix_put - use cvmix_kinds_and_types + ! CVMix is a fixed double-precision library (kind cvmix_r8). FESOM's WP may be real4 + ! (single precision); convert WP<->cvmix_r8 at the CVMix call boundary below so CVMix + ! always runs in double precision. In a double-precision FESOM build cvmix_r8==WP, so + ! the casts/temps are exact no-ops and the result is bit-identical. cvmix_r8 is provided + ! by the full use of cvmix_kinds_and_types below. + use cvmix_kinds_and_types !___________________________________________________________________________ ! module calls from FESOM @@ -316,7 +321,9 @@ subroutine init_cvmix_idemix(partit, mesh) !_______________________________________________________________________ ! initialise IDEMIX parameters - call init_idemix(idemix_tau_v,idemix_tau_h,idemix_gamma,idemix_jstar,idemix_mu0)! ,handle_old_vals)! ,idemix_userdef_constants) + call init_idemix(real(idemix_tau_v, cvmix_r8), real(idemix_tau_h, cvmix_r8), & + real(idemix_gamma, cvmix_r8), real(idemix_jstar, cvmix_r8), & + real(idemix_mu0, cvmix_r8))! ,handle_old_vals)! ,idemix_userdef_constants) end subroutine init_cvmix_idemix ! ! @@ -334,6 +341,16 @@ subroutine calc_cvmix_idemix(partit, mesh) real(kind=WP) :: grad_v0Eiw(2), deltaX1, deltaY1, deltaX2, deltaY2 real(kind=WP) :: tsum1, tsum2, tsum3, tsum4, tvol logical :: debug=.false. + ! double-precision (cvmix_r8) buffers for the cvmix_coeffs_idemix call boundary. + ! One per array argument; sliced identically to the original WP actual arguments so + ! the sequence association is unchanged. Exact no-op when WP==cvmix_r8. + real(cvmix_r8) :: dzw_r8(mesh%nl), dzt_r8(mesh%nl), iwe_new_r8(mesh%nl) + real(cvmix_r8) :: iwe_old_r8(mesh%nl), alpha_c_r8(mesh%nl) + real(cvmix_r8) :: KappaM_r8(mesh%nl), KappaH_r8(mesh%nl), Nsqr_r8(mesh%nl) + real(cvmix_r8) :: iwe_Ttot_r8(mesh%nl), iwe_Tdif_r8(mesh%nl), iwe_Thdi_r8(mesh%nl) + real(cvmix_r8) :: iwe_Tdis_r8(mesh%nl), iwe_Tsur_r8(mesh%nl), iwe_Tbot_r8(mesh%nl) + real(cvmix_r8) :: c0_r8(mesh%nl), v0_r8(mesh%nl) + real(cvmix_r8) :: cvmix_int_1_r8(mesh%nl), cvmix_int_2_r8(mesh%nl), cvmix_int_3_r8(mesh%nl) #include "../associate_part_def.h" #include "../associate_mesh_def.h" @@ -369,46 +386,72 @@ subroutine calc_cvmix_idemix(partit, mesh) !___________________________________________________________________ ! main call to calculate idemix iwe_old = iwe(:,elem) + !___________________________________________________________________ + ! copy WP input/inout arrays into cvmix_r8 temps before the call + ! (slices match the original actual arguments exactly) + dzw_r8(uln:nln) = real(helem(uln:nln, elem), cvmix_r8) + dzt_r8(uln:nln+1) = real(dz_trr(uln:nln+1), cvmix_r8) + iwe_old_r8(uln:nln+1) = real(iwe_old(uln:nln+1), cvmix_r8) + Nsqr_r8(uln:nln+1) = real(bvfreq2(uln:nln+1), cvmix_r8) + iwe_Thdi_r8(uln:nln+1) = real(iwe_Thdi(uln:nln+1,elem), cvmix_r8) + KappaM_r8(uln:nln+1) = real(iwe_Av(uln:nln+1, elem), cvmix_r8) ! inout + KappaH_r8(uln:nln+1) = real(iwe_Kv(uln:nln+1, elem), cvmix_r8) ! inout call cvmix_coeffs_idemix(& ! parameter - dzw = helem(uln:nln, elem), & - dzt = dz_trr(uln:nln+1), & + dzw = dzw_r8(uln:nln), & + dzt = dzt_r8(uln:nln+1), & ! nlev = nln, & nlev = nln-uln+1, & - max_nlev = nl-1, & - dtime = dt, & - coriolis = mesh%coriolis(elem), & - ! essentials - iwe_new = iwe(uln:nln+1,elem), & ! out - iwe_old = iwe_old(uln:nln+1), & ! in - forc_iw_surface = iwe_fsrf(elem), & ! in - forc_iw_bottom = iwe_fbot(elem), & ! in + max_nlev = nl-1, & + dtime = real(dt, cvmix_r8), & + coriolis = real(mesh%coriolis(elem), cvmix_r8),& + ! essentials + iwe_new = iwe_new_r8(uln:nln+1), & ! out + iwe_old = iwe_old_r8(uln:nln+1), & ! in + forc_iw_surface = real(iwe_fsrf(elem), cvmix_r8), & ! in + forc_iw_bottom = real(iwe_fbot(elem), cvmix_r8), & ! in ! FIXME: nils: better output IDEMIX Ri directly - alpha_c = iwe_alpha_c(uln:nln+1, elem), & ! out (for Ri IMIX) - ! only for Osborn shortcut + alpha_c = alpha_c_r8(uln:nln+1), & ! out (for Ri IMIX) + ! only for Osborn shortcut ! FIXME: nils: put this to cvmix_tke - KappaM_out = iwe_Av( uln:nln+1, elem), & ! out - KappaH_out = iwe_Kv( uln:nln+1, elem), & ! out - Nsqr = bvfreq2( uln:nln+1), & ! in + KappaM_out = KappaM_r8(uln:nln+1), & ! inout + KappaH_out = KappaH_r8(uln:nln+1), & ! inout + Nsqr = Nsqr_r8(uln:nln+1), & ! in ! diagnostics - iwe_Ttot = iwe_Ttot(uln:nln+1, elem), & - iwe_Tdif = iwe_Tdif(uln:nln+1, elem), & - iwe_Thdi = iwe_Thdi(uln:nln+1, elem), & - iwe_Tdis = iwe_Tdis(uln:nln+1, elem), & - iwe_Tsur = iwe_Tsur(uln:nln+1, elem), & - iwe_Tbot = iwe_Tbot(uln:nln+1, elem), & - c0 = iwe_c0( uln:nln+1, elem), & - v0 = iwe_v0( uln:nln+1, elem), & + iwe_Ttot = iwe_Ttot_r8(uln:nln+1), & + iwe_Tdif = iwe_Tdif_r8(uln:nln+1), & + iwe_Thdi = iwe_Thdi_r8(uln:nln+1), & + iwe_Tdis = iwe_Tdis_r8(uln:nln+1), & + iwe_Tsur = iwe_Tsur_r8(uln:nln+1), & + iwe_Tbot = iwe_Tbot_r8(uln:nln+1), & + c0 = c0_r8(uln:nln+1), & + v0 = v0_r8(uln:nln+1), & ! debugging debug = debug, & !i = i, & !j = j, & !tstep_count = tstep_count, & - cvmix_int_1 = cvmix_dummy_1(uln:nln+1, elem), & - cvmix_int_2 = cvmix_dummy_2(uln:nln+1, elem), & - cvmix_int_3 = cvmix_dummy_3(uln:nln+1, elem) & + cvmix_int_1 = cvmix_int_1_r8(uln:nln+1), & + cvmix_int_2 = cvmix_int_2_r8(uln:nln+1), & + cvmix_int_3 = cvmix_int_3_r8(uln:nln+1) & ) - + !___________________________________________________________________ + ! copy cvmix_r8 output/inout temps back into the WP FESOM arrays + iwe(uln:nln+1,elem) = real(iwe_new_r8(uln:nln+1), WP) + iwe_alpha_c(uln:nln+1,elem) = real(alpha_c_r8(uln:nln+1), WP) + iwe_Av(uln:nln+1, elem) = real(KappaM_r8(uln:nln+1), WP) + iwe_Kv(uln:nln+1, elem) = real(KappaH_r8(uln:nln+1), WP) + iwe_Ttot(uln:nln+1, elem) = real(iwe_Ttot_r8(uln:nln+1), WP) + iwe_Tdif(uln:nln+1, elem) = real(iwe_Tdif_r8(uln:nln+1), WP) + iwe_Tdis(uln:nln+1, elem) = real(iwe_Tdis_r8(uln:nln+1), WP) + iwe_Tsur(uln:nln+1, elem) = real(iwe_Tsur_r8(uln:nln+1), WP) + iwe_Tbot(uln:nln+1, elem) = real(iwe_Tbot_r8(uln:nln+1), WP) + iwe_c0(uln:nln+1, elem) = real(c0_r8(uln:nln+1), WP) + iwe_v0(uln:nln+1, elem) = real(v0_r8(uln:nln+1), WP) + cvmix_dummy_1(uln:nln+1,elem) = real(cvmix_int_1_r8(uln:nln+1), WP) + cvmix_dummy_2(uln:nln+1,elem) = real(cvmix_int_2_r8(uln:nln+1), WP) + cvmix_dummy_3(uln:nln+1,elem) = real(cvmix_int_3_r8(uln:nln+1), WP) + end do !-->do node = 1,node_size !_______________________________________________________________________ diff --git a/src/cvmix_driver/gen_modules_cvmix_kpp.F90 b/src/cvmix_driver/gen_modules_cvmix_kpp.F90 index 62d0ff1a8..15a7e8271 100644 --- a/src/cvmix_driver/gen_modules_cvmix_kpp.F90 +++ b/src/cvmix_driver/gen_modules_cvmix_kpp.F90 @@ -412,11 +412,17 @@ subroutine init_cvmix_kpp(partit, mesh) !_______________________________________________________________________ ! Initialise CVMIX ! call the cvmix subroutine to initialise all required namelists - call cvmix_init_kpp(Ri_crit = kpp_Rib_crit, & - minOBLdepth = kpp_minOBLdepth, & - minVtsqr = kpp_minVtsqr, & - vonKarman = kpp_vonKarman, & - surf_layer_ext = kpp_surf_layer_ext, & + ! CVMix is a fixed double-precision library (kind cvmix_r8). FESOM's WP may + ! be real4 (single precision); convert WP<->cvmix_r8 at the CVMix call + ! boundary below so CVMix always runs in double precision. In a double- + ! precision FESOM build cvmix_r8==WP, so the casts/temps are exact no-ops + ! and the result is bit-identical. (cvmix_r8 is provided by the + ! "use cvmix_kinds_and_types" above.) + call cvmix_init_kpp(Ri_crit = real(kpp_Rib_crit, cvmix_r8), & + minOBLdepth = real(kpp_minOBLdepth, cvmix_r8), & + minVtsqr = real(kpp_minVtsqr, cvmix_r8), & + vonKarman = real(kpp_vonKarman, cvmix_r8), & + surf_layer_ext = real(kpp_surf_layer_ext, cvmix_r8), & interp_type = kpp_interptype_ri, & interp_type2 = kpp_interptype_atobl, & lEkman = kpp_use_compEkman, & @@ -458,6 +464,23 @@ subroutine calc_cvmix_kpp(ice, dynamics, tracers, partit, mesh) real(kind=WP), dimension(:) , pointer :: a_ice real(kind=WP), dimension(:,:) , pointer :: temp, salt real(kind=WP), dimension(:,:,:), pointer :: UVnode + ! CVMix is a fixed double-precision library (kind cvmix_r8); FESOM's WP may be + ! real4. These cvmix_r8 buffers hold copy-in/copy-out of the WP actuals at the + ! CVMix call boundary so CVMix always runs in double precision. When WP==cvmix_r8 + ! (double-precision build) every copy is an exact no-op (bit-identical results). + ! All are sized mesh%nl (max levels); calls use index sections matching the WP arrays. + real(cvmix_r8) :: zbar_r8(mesh%nl), zcntr_r8(mesh%nl) ! interfaces / centers + real(cvmix_r8) :: uE_r8(mesh%nl), vE_r8(mesh%nl) ! Eulerian vel (StokesXi) + real(cvmix_r8) :: uS_r8(mesh%nl), vS_r8(mesh%nl) ! Stokes drift inout (StokesXi) + real(cvmix_r8) :: uSbar_r8(mesh%nl), vSbar_r8(mesh%nl) ! cell-avg Stokes inout (StokesXi) + real(cvmix_r8) :: obl_r8(mesh%nl), buoy_r8(mesh%nl) ! abs(Z) / surf buoy flux + real(cvmix_r8) :: xi_r8(mesh%nl), ws_r8(mesh%nl) ! Stokes xi / turb scale w_s + real(cvmix_r8) :: dbuoy_r8(mesh%nl), dvsqr_r8(mesh%nl) ! bulk Ri inputs + real(cvmix_r8) :: nsqr_r8(mesh%nl), bulkRi_r8(mesh%nl) ! N^2 / bulk Ri result + real(cvmix_r8) :: ribulk_r8(mesh%nl), surfbuoy_r8(mesh%nl) ! OBL-depth inputs + real(cvmix_r8) :: md_r8(mesh%nl), td_r8(mesh%nl), sd_r8(mesh%nl) ! coeffs_kpp M/T/S diff (in+out) + real(cvmix_r8) :: tnl_r8(mesh%nl), snl_r8(mesh%nl) ! coeffs_kpp nonlocal T/S (in+out) + real(cvmix_r8) :: stokesXi_r8, obldep_r8, kobl_r8 ! scalar outputs #include "../associate_part_def.h" #include "../associate_mesh_def.h" #include "../associate_part_ass.h" @@ -603,29 +626,45 @@ subroutine calc_cvmix_kpp(ice, dynamics, tracers, partit, mesh) uS_sld_m, vS_sld_m) + ! --- CVMix boundary: copy WP actuals into cvmix_r8 buffers (no-op when WP==cvmix_r8) + zbar_r8(nun:nln+1) = real(zbar_3d_n(nun:nln+1,node), cvmix_r8) + zcntr_r8(nun:nln) = real(Z_3d_n( nun:nln, node), cvmix_r8) + uE_r8(nun:nln) = real(UVnode(1, nun:nln, node), cvmix_r8) + vE_r8(nun:nln) = real(UVnode(2, nun:nln, node), cvmix_r8) + uS_r8(nun:nln) = real(kpp_uS_t( nun:nln), cvmix_r8) ! inout + vS_r8(nun:nln) = real(kpp_VS_t( nun:nln), cvmix_r8) ! inout + uSbar_r8(nun:nln) = real(kpp_uS_m( nun:nln), cvmix_r8) ! inout + vSbar_r8(nun:nln) = real(kpp_VS_m( nun:nln), cvmix_r8) ! inout + stokesXi_r8 = real(kpp_stokesXi_z(nz), cvmix_r8) ! inout (out) call cvmix_kpp_compute_StokesXi (& - zbar_3d_n(nun:nln+1,node), & ! full depth levels - Z_3d_n(nun:nln,node), & ! mid depth levels + zbar_r8(nun:nln+1), & ! full depth levels + zcntr_r8(nun:nln), & ! mid depth levels nzsfc, & ! cell index of Surface Layer Depth - sldepth, & ! surface layer depth > 0 - kpp_buoyflx_nl(nz), & ! surfce buoyancy flux (m2/s3) total - ustar, & ! turbulent friction velocity at surface (m/s), - mesh%coriolis_node(node), & ! Coriolis parameter (1/s) dim=1 - UVnode(1,nun:nln,node), & ! zonal Eulerian mean horizontal velocity components - UVnode(2,nun:nln,node), & ! merid Eulerian mean horizontal velocity components + real(sldepth, cvmix_r8), & ! surface layer depth > 0 + real(kpp_buoyflx_nl(nz), cvmix_r8), & ! surfce buoyancy flux (m2/s3) total + real(ustar, cvmix_r8), & ! turbulent friction velocity at surface (m/s), + real(mesh%coriolis_node(node), cvmix_r8), & ! Coriolis parameter (1/s) dim=1 + uE_r8(nun:nln), & ! zonal Eulerian mean horizontal velocity components + vE_r8(nun:nln), & ! merid Eulerian mean horizontal velocity components !___stokes velocities___________________________________________ - kpp_uS_t(nun:nln), & ! zonal Surface Stokes drift velocity (Wave-induced drift) - kpp_VS_t(nun:nln), & ! merid Surface Stokes drift velocity (Wave-induced drift) - kpp_uS_m(nun:nln), & ! - kpp_VS_m(nun:nln), & ! - uS_sld_t, & ! - vS_sld_t, & ! - uS_sld_m, & ! - vS_sld_m, & ! + uS_r8(nun:nln), & ! zonal Surface Stokes drift velocity (Wave-induced drift) + vS_r8(nun:nln), & ! merid Surface Stokes drift velocity (Wave-induced drift) + uSbar_r8(nun:nln), & ! + vSbar_r8(nun:nln), & ! + real(uS_sld_t, cvmix_r8), & ! + real(vS_sld_t, cvmix_r8), & ! + real(uS_sld_m, cvmix_r8), & ! + real(vS_sld_m, cvmix_r8), & ! !___outputs_____________________________________________________ - kpp_stokesXi_z(nz) & ! (out) Stokes Similartiy parameter + stokesXi_r8 & ! (out) Stokes Similartiy parameter ) - + ! --- copy cvmix_r8 results back to WP (StokesXi restores the inout drift arrays) + kpp_uS_t(nun:nln) = real(uS_r8(nun:nln), WP) + kpp_VS_t(nun:nln) = real(vS_r8(nun:nln), WP) + kpp_uS_m(nun:nln) = real(uSbar_r8(nun:nln), WP) + kpp_VS_m(nun:nln) = real(vSbar_r8(nun:nln), WP) + kpp_stokesXi_z(nz) = real(stokesXi_r8, WP) + kpp_stokesVt_z(nz) = 0.0! kpp_stokesXi_z(nz) ! kpp_stokesXi_z: It represents the strength of Langmuir turbulence — @@ -793,14 +832,20 @@ subroutine calc_cvmix_kpp(ice, dynamics, tracers, partit, mesh) ! --> w_s is re-computed in call cvmix_coeffs_kpp_low under consideration ! StokesXi when CVmix_kpp_params_in%lStokesMOST = .True. ! --> sigma_coord must be here 1.0!!! --> see MOM6 --> parameterizations/vertical/MOM_CVMix_KPP.F90:1343 + ! --- CVMix boundary: copy WP actuals into cvmix_r8 buffers (no-op when WP==cvmix_r8) + obl_r8(nun:nln) = real(abs(Z_3d_n(nun:nln,node)), cvmix_r8) + buoy_r8(nun:nln) = real(kpp_buoyflx_nl(nun:nln), cvmix_r8) + xi_r8(nun:nln) = real(kpp_stokesVt_z(nun:nln), cvmix_r8) + ws_r8(nun:nln) = real(kpp_ws_cntr(nun:nln), cvmix_r8) ! inout call cvmix_kpp_compute_turbulent_scales( & - sigma_coord = 1.0_WP , & ! (in) sigma: Normalized surface layer depth - OBL_depth = abs(Z_3d_n(nun:nln,node)), & ! (in) Assume OBL depth (m) = mid-depth level - surf_buoy_force = kpp_buoyflx_nl(nun:nln), & ! (in) surfce buoyancy flux (m2/s3) consider sw_pene - surf_fric_vel = ustar, & ! (in) turbulent friction velocity at surface (m/s) - xi = kpp_stokesVt_z(nun:nln), & ! (in) Stokes parameter xi= Ps/(PU+PS+PB) - w_s = kpp_ws_cntr(nun:nln) & ! (out) Turbulent velocity scale profile (m/s) for skalars - ) + sigma_coord = 1.0_cvmix_r8 , & ! (in) sigma: Normalized surface layer depth + OBL_depth = obl_r8(nun:nln), & ! (in) Assume OBL depth (m) = mid-depth level + surf_buoy_force = buoy_r8(nun:nln), & ! (in) surfce buoyancy flux (m2/s3) consider sw_pene + surf_fric_vel = real(ustar, cvmix_r8), & ! (in) turbulent friction velocity at surface (m/s) + xi = xi_r8(nun:nln), & ! (in) Stokes parameter xi= Ps/(PU+PS+PB) + w_s = ws_r8(nun:nln) & ! (out) Turbulent velocity scale profile (m/s) for skalars + ) + kpp_ws_cntr(nun:nln) = real(ws_r8(nun:nln), WP) ! --> need w_s to compute cvmix_kpp_compute_bulk_Richardson(...) @@ -818,7 +863,7 @@ subroutine calc_cvmix_kpp(ice, dynamics, tracers, partit, mesh) ! computes the surface layer averaged Stokes drift, given ! the 10-meter wind (m/s) and the boundary layer depth (m). - uv_SLmean = cvmix_kpp_ustokes_SL_model(wind_norm, oblguess, CVmix_params_in) + uv_SLmean = real(cvmix_kpp_ustokes_SL_model(real(wind_norm, cvmix_r8), real(oblguess, cvmix_r8), CVmix_params_in), WP) ! Copute Langmuir enhance,ent factor & SL langmuir number call compute_Efactor( & @@ -845,17 +890,26 @@ subroutine calc_cvmix_kpp(ice, dynamics, tracers, partit, mesh) ! --> v_t = turbulent velocity shear, needs to be parameterized ! --> v_t = z * ws * N * v_tc ! |-> v_tc = Cv * sqrt(0.2/Cs/epsilon) / kpp_vonKarman^2 / kpp_Rib_crit - kpp_bulkRi(nun:nln+1) = cvmix_kpp_compute_bulk_Richardson( & - zt_cntr = Z_3d_n( nun:nln ,node), & ! (in) Depth of cell center (m) - delta_buoy_cntr = kpp_dbsurf( nun:nln), & ! (in) Bulk buoyancy difference, Br-B(z) (1/s) - delta_Vsqr_cntr = kpp_dvsurf2(nun:nln), & ! (in) Square of resolved velocity difference (m2/s2) - ws_cntr = kpp_ws_cntr(nun:nln), & ! (in) Turbulent velocity scale profile (m/s) - Nsqr_iface = bvfreq( nun:nln+1,node), & ! (in) Buoyancy frequency (1/s) - bfsfc = kpp_buoyflx_nl(nun:nln), & ! (in) surface buoyancy flux (units: m^2/s^3) - ustar = ustar, & ! (in) friction velocity (units: m/s) - EFactor = kpp_EFactor(node), & ! (in) Langmuir enhancement factor for entrainment (units: none) - LaSL = kpp_LaSL(node) & ! (in) surface layer averaged Langmuir number (units: none) - ) + ! --- CVMix boundary: copy WP actuals into cvmix_r8 buffers (no-op when WP==cvmix_r8) + zcntr_r8(nun:nln) = real(Z_3d_n( nun:nln ,node), cvmix_r8) + dbuoy_r8(nun:nln) = real(kpp_dbsurf( nun:nln), cvmix_r8) + dvsqr_r8(nun:nln) = real(kpp_dvsurf2(nun:nln), cvmix_r8) + ws_r8(nun:nln) = real(kpp_ws_cntr(nun:nln), cvmix_r8) + nsqr_r8(nun:nln+1) = real(bvfreq( nun:nln+1,node), cvmix_r8) + buoy_r8(nun:nln) = real(kpp_buoyflx_nl(nun:nln), cvmix_r8) + ! NOTE: LHS section nun:nln+1 kept identical to the original to preserve behavior + bulkRi_r8(nun:nln+1) = cvmix_kpp_compute_bulk_Richardson( & + zt_cntr = zcntr_r8(nun:nln), & ! (in) Depth of cell center (m) + delta_buoy_cntr = dbuoy_r8(nun:nln), & ! (in) Bulk buoyancy difference, Br-B(z) (1/s) + delta_Vsqr_cntr = dvsqr_r8(nun:nln), & ! (in) Square of resolved velocity difference (m2/s2) + ws_cntr = ws_r8(nun:nln), & ! (in) Turbulent velocity scale profile (m/s) + Nsqr_iface = nsqr_r8(nun:nln+1), & ! (in) Buoyancy frequency (1/s) + bfsfc = buoy_r8(nun:nln), & ! (in) surface buoyancy flux (units: m^2/s^3) + ustar = real(ustar, cvmix_r8), & ! (in) friction velocity (units: m/s) + EFactor = real(kpp_EFactor(node), cvmix_r8), & ! (in) Langmuir enhancement factor for entrainment (units: none) + LaSL = real(kpp_LaSL(node), cvmix_r8) & ! (in) surface layer averaged Langmuir number (units: none) + ) + kpp_bulkRi(nun:nln+1) = real(bulkRi_r8(nun:nln+1), WP) !___________________________________________________________________ @@ -873,17 +927,25 @@ subroutine calc_cvmix_kpp(ice, dynamics, tracers, partit, mesh) ! everywhere and eradicates possible obldepths from bulkRi number ! --> thats why kpp_use_compEkman = .False. ! --> check if its still the case + ! --- CVMix boundary: copy WP actuals into cvmix_r8 buffers (no-op when WP==cvmix_r8) + ribulk_r8(nun:nln+1) = real(kpp_bulkRi(nun:nln+1), cvmix_r8) + zbar_r8(nun:nln+1) = real(zbar_3d_n(nun:nln+1,node), cvmix_r8) + zcntr_r8(nun:nln) = real(Z_3d_n( nun:nln, node), cvmix_r8) + surfbuoy_r8(nun:nln) = real(aux_surfbuoyflx_nl(nun:nln), cvmix_r8) + xi_r8(nun:nln) = real(kpp_stokesXi_z(nun:nln), cvmix_r8) call cvmix_kpp_compute_OBL_depth( & - Ri_bulk = kpp_bulkRi(nun:nln+1), & ! (in) Bulk Richardson number dim=(ke+1) - zw_iface = zbar_3d_n( nun:nln+1,node), & ! (in) Height of interfaces (m) dim=(ke+1) - OBL_depth = kpp_obldepth(node), & ! (out) OBL depth (m) dim=1 - kOBL_depth = kpp_nzobldepth(node), & ! (out) level (+fraction) of OBL extent dim=1 - zt_cntr = Z_3d_n( nun:nln ,node), & ! (in) Height of cell centers (m) dim=(ke) - surf_fric = ustar, & ! (in) Turbulent friction velocity at surface (m/s) dim=1 - surf_buoy = aux_surfbuoyflx_nl(nun:nln),& ! (in) Buoyancy flux at surface (m2/s3) dim=1 - Coriolis = mesh%coriolis_node(node), & ! (in) Coriolis parameter (1/s) dim=1 - Xi = kpp_stokesXi_z(nun:nln) & - ) + Ri_bulk = ribulk_r8(nun:nln+1), & ! (in) Bulk Richardson number dim=(ke+1) + zw_iface = zbar_r8(nun:nln+1), & ! (in) Height of interfaces (m) dim=(ke+1) + OBL_depth = obldep_r8, & ! (out) OBL depth (m) dim=1 + kOBL_depth = kobl_r8, & ! (out) level (+fraction) of OBL extent dim=1 + zt_cntr = zcntr_r8(nun:nln), & ! (in) Height of cell centers (m) dim=(ke) + surf_fric = real(ustar, cvmix_r8), & ! (in) Turbulent friction velocity at surface (m/s) dim=1 + surf_buoy = surfbuoy_r8(nun:nln), & ! (in) Buoyancy flux at surface (m2/s3) dim=1 + Coriolis = real(mesh%coriolis_node(node), cvmix_r8), & ! (in) Coriolis parameter (1/s) dim=1 + Xi = xi_r8(nun:nln) & + ) + kpp_obldepth(node) = real(obldep_r8, WP) + kpp_nzobldepth(node) = real(kobl_r8, WP) kpp_nzobldepth(node) = kpp_nzobldepth(node) + nun - 1 !___safty switches for kpp_nzobldepth_______________________________ @@ -951,32 +1013,48 @@ subroutine calc_cvmix_kpp(ice, dynamics, tracers, partit, mesh) ! Cray compiler workaround: store subscript as scalar integer nz_buoyflx_idx = kpp_nzobldepth(node) + ! --- CVMix boundary: copy WP actuals into cvmix_r8 buffers (no-op when WP==cvmix_r8) + zbar_r8(nun:nln+1) = real(zbar_3d_n(nun:nln+1,node), cvmix_r8) + zcntr_r8(nun:nln) = real(Z_3d_n( nun:nln, node), cvmix_r8) + uE_r8(nun:nln) = real(UVnode(1, nun:nln, node), cvmix_r8) + vE_r8(nun:nln) = real(UVnode(2, nun:nln, node), cvmix_r8) + uS_r8(nun:nln) = real(kpp_uS_t( nun:nln), cvmix_r8) ! inout + vS_r8(nun:nln) = real(kpp_VS_t( nun:nln), cvmix_r8) ! inout + uSbar_r8(nun:nln) = real(kpp_uS_m( nun:nln), cvmix_r8) ! inout + vSbar_r8(nun:nln) = real(kpp_VS_m( nun:nln), cvmix_r8) ! inout + stokesXi_r8 = real(kpp_stokesXi(node), cvmix_r8) ! inout (out) call cvmix_kpp_compute_StokesXi (& - zbar_3d_n(nun:nln+1,node), & ! full depth levels - Z_3d_n(nun:nln,node), & ! mid depth levels + zbar_r8(nun:nln+1), & ! full depth levels + zcntr_r8(nun:nln), & ! mid depth levels nzsfc, & ! cell index of Surface Layer Depth - sldepth, & ! surface layer depth > 0 - kpp_buoyflx_nl(nz_buoyflx_idx), & ! surfce buoyancy flux (m2/s3) total - ustar, & ! turbulent friction velocity at surface (m/s), - mesh%coriolis_node(node), & ! Coriolis parameter (1/s) dim=1 - UVnode(1,nun:nln,node), & ! zonal Eulerian mean horizontal velocity components - UVnode(2,nun:nln,node), & ! merid Eulerian mean horizontal velocity components + real(sldepth, cvmix_r8), & ! surface layer depth > 0 + real(kpp_buoyflx_nl(nz_buoyflx_idx), cvmix_r8), & ! surfce buoyancy flux (m2/s3) total + real(ustar, cvmix_r8), & ! turbulent friction velocity at surface (m/s), + real(mesh%coriolis_node(node), cvmix_r8), & ! Coriolis parameter (1/s) dim=1 + uE_r8(nun:nln), & ! zonal Eulerian mean horizontal velocity components + vE_r8(nun:nln), & ! merid Eulerian mean horizontal velocity components !___stokes velocities___________________________________________ - kpp_uS_t(nun:nln), & ! zonal Surface Stokes drift velocity (Wave-induced drift) - kpp_VS_t(nun:nln), & ! merid Surface Stokes drift velocity (Wave-induced drift) - kpp_uS_m(nun:nln), & ! - kpp_VS_m(nun:nln), & ! - uS_sld_t, & ! - vS_sld_t, & ! - uS_sld_m, & ! - vS_sld_m, & ! + uS_r8(nun:nln), & ! zonal Surface Stokes drift velocity (Wave-induced drift) + vS_r8(nun:nln), & ! merid Surface Stokes drift velocity (Wave-induced drift) + uSbar_r8(nun:nln), & ! + vSbar_r8(nun:nln), & ! + real(uS_sld_t, cvmix_r8), & ! + real(vS_sld_t, cvmix_r8), & ! + real(uS_sld_m, cvmix_r8), & ! + real(vS_sld_m, cvmix_r8), & ! !___outputs_____________________________________________________ - kpp_stokesXi(node) ) ! (out) Stokes Similartiy parameter + stokesXi_r8 ) ! (out) Stokes Similartiy parameter + ! --- copy cvmix_r8 results back to WP (StokesXi restores the inout drift arrays) + kpp_uS_t(nun:nln) = real(uS_r8(nun:nln), WP) + kpp_VS_t(nun:nln) = real(vS_r8(nun:nln), WP) + kpp_uS_m(nun:nln) = real(uSbar_r8(nun:nln), WP) + kpp_VS_m(nun:nln) = real(vSbar_r8(nun:nln), WP) + kpp_stokesXi(node) = real(stokesXi_r8, WP) end if ! --> if (kpp_use_StokesMOST) then ! computes the surface layer averaged Stokes drift, given ! the 10-meter wind (m/s) and the boundary layer depth (m). - uv_SLmean = cvmix_kpp_ustokes_SL_model(wind_norm, kpp_obldepth(node), CVmix_params_in) + uv_SLmean = real(cvmix_kpp_ustokes_SL_model(real(wind_norm, cvmix_r8), real(kpp_obldepth(node), cvmix_r8), CVmix_params_in), WP) ! Copute Langmuir enhance,ent factor & SL langmuir number call compute_Efactor( & @@ -1053,26 +1131,42 @@ subroutine calc_cvmix_kpp(ice, dynamics, tracers, partit, mesh) !___________________________________________________________________ ! 11) compute the turbulent diffusion coefficients !!PS if (flag_debug .and. mype==0) print *, achar(27)//'[35m'//' --> calc kpp coeff'//achar(27)//'[0m' + ! --- CVMix boundary: copy WP actuals into cvmix_r8 buffers (no-op when WP==cvmix_r8) + ! md/td/sd_r8 serve BOTH the *_out (inout) and old_* (in) args (same actual as the + ! original kpp_oblmixc aliasing); seed them with the current values, then copy out. + md_r8(:) = real(kpp_oblmixc(:,node,1), cvmix_r8) ! Mdiff_out (inout) & old_Mdiff (in) + td_r8(:) = real(kpp_oblmixc(:,node,2), cvmix_r8) ! Tdiff_out (inout) & old_Tdiff (in) + sd_r8(:) = real(kpp_oblmixc(:,node,3), cvmix_r8) ! Sdiff_out (inout) & old_Sdiff (in) + zbar_r8(:) = real(zbar_3d_n(:,node), cvmix_r8) ! zw (in) + zcntr_r8(:) = real(Z_3d_n(:,node), cvmix_r8) ! zt (in) + tnl_r8(:) = real(kpp_nonlcltranspT(:,node), cvmix_r8) ! Tnonlocal (inout) + snl_r8(:) = real(kpp_nonlcltranspS(:,node), cvmix_r8) ! Snonlocal (inout) call cvmix_coeffs_kpp( & - Mdiff_out = kpp_oblmixc(:,node,1), & ! (inout) new_Mdiff: Total viscosity (m2/s) - Tdiff_out = kpp_oblmixc(:,node,2), & ! (inout) new_Tdiff: Total heat diffusivity (m2/s) - Sdiff_out = kpp_oblmixc(:,node,3), & ! (inout) new_Sdiff: Total salt diffusivity (m2/s) - zw = zbar_3d_n(:,node), & ! (in) zw_iface: Height of interfaces (m) - zt = Z_3d_n(:,node), & ! (in) zt_cntr: Height of level centers (m) - old_Mdiff = kpp_oblmixc(:,node,1), & ! (in) MDiff_iface: Original viscosity (m2/s) - old_Tdiff = kpp_oblmixc(:,node,2), & ! (in) Tdiff_iface: Original heat diffusivity (m2/s) - old_Sdiff = kpp_oblmixc(:,node,3), & ! (in) Sdiff_iface: Original salt diffusivity (m2/s) - OBL_depth = kpp_obldepth(node), & ! (in) BoundaryLayerDepth: OBL depth (m) - kOBL_depth= kpp_nzobldepth(node), & ! (in) kOBL_depth: level (+fraction) of OBL extent - Tnonlocal = kpp_nonlcltranspT(:,node), & ! (out) kpp_Tnonlocal_iface: Non-local heat transport (non-dimensional) - Snonlocal = kpp_nonlcltranspS(:,node), & ! (out) kkp_Snonlocal_iface: Non-local salt transport (non-dimensional) - surf_fric = ustar, & ! (in) SurfaceFriction:Turbulent friction velocity at surface (m/s) - surf_buoy = aux_surfbuoyflx_nl(1), & ! (in) SurfaceBuoynacyForcing: Buoyancy flux at surface (m2/s3) + Mdiff_out = md_r8, & ! (inout) new_Mdiff: Total viscosity (m2/s) + Tdiff_out = td_r8, & ! (inout) new_Tdiff: Total heat diffusivity (m2/s) + Sdiff_out = sd_r8, & ! (inout) new_Sdiff: Total salt diffusivity (m2/s) + zw = zbar_r8, & ! (in) zw_iface: Height of interfaces (m) + zt = zcntr_r8, & ! (in) zt_cntr: Height of level centers (m) + old_Mdiff = md_r8, & ! (in) MDiff_iface: Original viscosity (m2/s) + old_Tdiff = td_r8, & ! (in) Tdiff_iface: Original heat diffusivity (m2/s) + old_Sdiff = sd_r8, & ! (in) Sdiff_iface: Original salt diffusivity (m2/s) + OBL_depth = real(kpp_obldepth(node), cvmix_r8), & ! (in) BoundaryLayerDepth: OBL depth (m) + kOBL_depth= real(kpp_nzobldepth(node), cvmix_r8), & ! (in) kOBL_depth: level (+fraction) of OBL extent + Tnonlocal = tnl_r8, & ! (out) kpp_Tnonlocal_iface: Non-local heat transport (non-dimensional) + Snonlocal = snl_r8, & ! (out) kkp_Snonlocal_iface: Non-local salt transport (non-dimensional) + surf_fric = real(ustar, cvmix_r8), & ! (in) SurfaceFriction:Turbulent friction velocity at surface (m/s) + surf_buoy = real(aux_surfbuoyflx_nl(1), cvmix_r8), & ! (in) SurfaceBuoynacyForcing: Buoyancy flux at surface (m2/s3) nlev = nlevels_nod2D(node)-1, & ! (in) nlev: Number of levels to compute coeffs for max_nlev = nl-1, & ! (in) max_lev: maximum vertical levels - Langmuir_EFactor = kpp_EFactor(node), & ! (in) Langmuir enhancement factor - StokesXI = kpp_stokesXI(node) & ! + Langmuir_EFactor = real(kpp_EFactor(node), cvmix_r8), & ! (in) Langmuir enhancement factor + StokesXI = real(kpp_stokesXI(node), cvmix_r8) & ! ) + ! --- copy cvmix_r8 results back to WP + kpp_oblmixc(:,node,1) = real(md_r8, WP) + kpp_oblmixc(:,node,2) = real(td_r8, WP) + kpp_oblmixc(:,node,3) = real(sd_r8, WP) + kpp_nonlcltranspT(:,node) = real(tnl_r8, WP) + kpp_nonlcltranspS(:,node) = real(snl_r8, WP) !___________________________________________________________________ diff --git a/src/cvmix_driver/gen_modules_cvmix_pp.F90 b/src/cvmix_driver/gen_modules_cvmix_pp.F90 index 478b6b7f6..7878f4884 100644 --- a/src/cvmix_driver/gen_modules_cvmix_pp.F90 +++ b/src/cvmix_driver/gen_modules_cvmix_pp.F90 @@ -19,7 +19,12 @@ module g_cvmix_pp !___________________________________________________________________________ ! module calls from cvmix library use cvmix_shear, only: cvmix_init_shear, cvmix_coeffs_shear - + ! CVMix is a fixed double-precision library (kind cvmix_r8). FESOM's WP may be real4 + ! (single precision); convert WP<->cvmix_r8 at the CVMix call boundary below so CVMix + ! always runs in double precision. In a double-precision FESOM build cvmix_r8==WP, so + ! the casts/temps are exact no-ops and the result is bit-identical. + use cvmix_kinds_and_types, only: cvmix_r8 + !___________________________________________________________________________ ! module calls from FESOM use g_config @@ -142,19 +147,19 @@ subroutine init_cvmix_pp(partit, mesh) ! already diffusive model even more diffusive --> it was done in ! FESOM1.4 like this. In this case set pp_Avbckg and pp_Kvbckg by ! hand - call cvmix_init_shear(mix_scheme = 'PP', & - PP_nu_zero = pp_Av0, & - PP_alpha = pp_alpha, & - PP_exp = pp_exp, & - PP_nu_b = 0.0_WP, & - PP_kappa_b = 0.0_WP) + call cvmix_init_shear(mix_scheme = 'PP', & + PP_nu_zero = real(pp_Av0, cvmix_r8), & + PP_alpha = real(pp_alpha, cvmix_r8), & + PP_exp = real(pp_exp, cvmix_r8), & + PP_nu_b = 0.0_cvmix_r8, & + PP_kappa_b = 0.0_cvmix_r8) else - call cvmix_init_shear(mix_scheme = 'PP', & - PP_nu_zero = pp_Av0, & - PP_alpha = pp_alpha, & - PP_exp = pp_exp, & - PP_nu_b = pp_Avbckg, & - PP_kappa_b = pp_Kvbckg) + call cvmix_init_shear(mix_scheme = 'PP', & + PP_nu_zero = real(pp_Av0, cvmix_r8), & + PP_alpha = real(pp_alpha, cvmix_r8), & + PP_exp = real(pp_exp, cvmix_r8), & + PP_nu_b = real(pp_Avbckg, cvmix_r8), & + PP_kappa_b = real(pp_Kvbckg, cvmix_r8)) end if end subroutine init_cvmix_pp ! @@ -171,6 +176,8 @@ subroutine calc_cvmix_pp(dynamics, partit, mesh) type(t_dyn), intent(inout), target :: dynamics integer :: node, elem, nz, nln, nun, elnodes(3), windnl=2, node_size real(kind=WP) :: vshear2, dz2, Kvb + ! double-precision (cvmix_r8) buffers for the CVMix call boundary + real(cvmix_r8) :: av_r8(mesh%nl), kv_r8(mesh%nl), rich_r8(mesh%nl) real(kind=WP), dimension(:,:,:), pointer :: UVnode #include "../associate_part_def.h" #include "../associate_mesh_def.h" @@ -212,11 +219,14 @@ subroutine calc_cvmix_pp(dynamics, partit, mesh) !___________________________________________________________________ ! use cvmix library function - call cvmix_coeffs_shear(Mdiff_out = pp_Av(:,node), & - Tdiff_out = pp_Kv(:,node), & - RICH = pp_richardnmb(:,node), & + rich_r8(:) = real(pp_richardnmb(:,node), cvmix_r8) + call cvmix_coeffs_shear(Mdiff_out = av_r8, & + Tdiff_out = kv_r8, & + RICH = rich_r8, & nlev = nln, & max_nlev = nl-1) + pp_Av(:,node) = real(av_r8, WP) + pp_Kv(:,node) = real(kv_r8, WP) !___________________________________________________________________ ! In the fesom flavour of PP mixing is the term from the background diff --git a/src/cvmix_driver/gen_modules_cvmix_tidal.F90 b/src/cvmix_driver/gen_modules_cvmix_tidal.F90 index c61a841c9..2628340c4 100644 --- a/src/cvmix_driver/gen_modules_cvmix_tidal.F90 +++ b/src/cvmix_driver/gen_modules_cvmix_tidal.F90 @@ -8,7 +8,7 @@ module g_cvmix_tidal !___________________________________________________________________________ ! module calls from cvmix library use cvmix_tidal, only : cvmix_init_tidal, cvmix_coeffs_tidal, cvmix_compute_Simmons_invariant - use cvmix_kinds_and_types, only: cvmix_data_type, cvmix_global_params_type + use cvmix_kinds_and_types, only: cvmix_data_type, cvmix_global_params_type, cvmix_r8 !___________________________________________________________________________ ! module calls from FESOM @@ -204,12 +204,12 @@ subroutine init_cvmix_tidal(partit, mesh) !_______________________________________________________________________ ! initialise TIDAL parameters - call cvmix_init_tidal(mix_scheme = tidal_mixscheme, & - efficiency = tidal_efficiency, & - vertical_decay_scale = tidal_vert_decayscale, & - max_coefficient = tidal_max_coeff, & - local_mixing_frac = tidal_lcl_mixfrac, & - depth_cutoff = tidal_depth_cutoff) + call cvmix_init_tidal(mix_scheme = tidal_mixscheme, & + efficiency = real(tidal_efficiency, cvmix_r8), & + vertical_decay_scale = real(tidal_vert_decayscale, cvmix_r8), & + max_coefficient = real(tidal_max_coeff, cvmix_r8), & + local_mixing_frac = real(tidal_lcl_mixfrac, cvmix_r8), & + depth_cutoff = real(tidal_depth_cutoff, cvmix_r8)) end subroutine init_cvmix_tidal ! @@ -226,6 +226,11 @@ subroutine calc_cvmix_tidal(partit, mesh) real(kind=WP) :: simmonscoeff, vertdep(mesh%nl) real(kind=WP) :: zbar_e(mesh%nl), Z_e(mesh%nl-1), bvfreq2(mesh%nl) real(kind=WP) :: tsum1, tvol + ! cvmix_r8 temporaries so CVMix always runs in double precision (no-op when WP==cvmix_r8) + real(cvmix_r8) :: simmonscoeff_r8 + real(cvmix_r8) :: vertdep_r8(mesh%nl), bvfreq2_r8(mesh%nl) + real(cvmix_r8) :: zbar_e_r8(mesh%nl), Z_e_r8(mesh%nl-1) + real(cvmix_r8) :: av_r8(mesh%nl), kv_r8(mesh%nl) #include "../associate_part_def.h" #include "../associate_mesh_def.h" #include "../associate_part_ass.h" @@ -263,28 +268,41 @@ subroutine calc_cvmix_tidal(partit, mesh) !___________________________________________________________________ ! Compute the time-invariant portion of the tidal mixing coefficient ! using the Simmons et al.(2004) scheme. + ! copy WP -> cvmix_r8 for input arrays before the CVMix call + zbar_e_r8(uln:nln+1) = real(zbar_e(uln:nln+1), cvmix_r8) + Z_e_r8(uln:nln) = real(Z_e(uln:nln), cvmix_r8) + vertdep_r8(uln:nln) = real(vertdep(uln:nln), cvmix_r8) call cvmix_compute_Simmons_invariant( & nlev = nln-uln+1 , & - energy_flux = tidal_fbot(elem), & !in W m-2 - rho = density_0, & - SimmonsCoeff = simmonscoeff, & - VertDep = vertdep(uln:nln), & ! vertical deposition function - zw = zbar_e(uln:nln+1), & - zt = Z_e(uln:nln)) + energy_flux = real(tidal_fbot(elem), cvmix_r8), & !in W m-2 + rho = real(density_0, cvmix_r8), & + SimmonsCoeff = simmonscoeff_r8, & + VertDep = vertdep_r8(uln:nln), & ! vertical deposition function + zw = zbar_e_r8(uln:nln+1), & + zt = Z_e_r8(uln:nln)) + ! copy cvmix_r8 -> WP for output args after the CVMix call + simmonscoeff = real(simmonscoeff_r8, WP) + vertdep(uln:nln) = real(vertdep_r8(uln:nln), WP) !___________________________________________________________________ ! Computes vertical diffusion coefficients for tidal mixing ! parameterizations. + ! copy WP -> cvmix_r8 for input arrays before the CVMix call + bvfreq2_r8(uln:nln) = real(bvfreq2(uln:nln), cvmix_r8) + vertdep_r8(uln:nln) = real(vertdep(uln:nln), cvmix_r8) call cvmix_coeffs_tidal( & - Mdiff_out = tidal_Av(uln:nln,elem), & - Tdiff_out = tidal_Kv(uln:nln,elem), & - Nsqr = bvfreq2(uln:nln), & !FIXME: limit to N2 > 10^-8 ? as in Simmons et al. - OceanDepth = -zbar_e_bot(elem), & !FIXME: neglecting free surface contribution - SimmonsCoeff = simmonscoeff, & - vert_dep = vertdep(uln:nln), & + Mdiff_out = av_r8(uln:nln), & + Tdiff_out = kv_r8(uln:nln), & + Nsqr = bvfreq2_r8(uln:nln), & !FIXME: limit to N2 > 10^-8 ? as in Simmons et al. + OceanDepth = real(-zbar_e_bot(elem), cvmix_r8), & !FIXME: neglecting free surface contribution + SimmonsCoeff = real(simmonscoeff, cvmix_r8), & + vert_dep = vertdep_r8(uln:nln), & nlev = nln-uln+1, & max_nlev = nl-1, & CVmix_params = CVmix_tidal_params) ! FIXME: Simmons et al. use Prandtl=10.0 (atm its 1.0) + ! copy cvmix_r8 -> WP for output arrays after the CVMix call + tidal_Av(uln:nln,elem) = real(av_r8(uln:nln), WP) + tidal_Kv(uln:nln,elem) = real(kv_r8(uln:nln), WP) !!PS if (mype==0) then !!PS write(*,*) 'SimmonsCoeff = ',simmonscoeff diff --git a/src/cvmix_driver/gen_modules_cvmix_tke.F90 b/src/cvmix_driver/gen_modules_cvmix_tke.F90 index 44bdf6866..3c439a901 100644 --- a/src/cvmix_driver/gen_modules_cvmix_tke.F90 +++ b/src/cvmix_driver/gen_modules_cvmix_tke.F90 @@ -18,6 +18,8 @@ module g_cvmix_tke ! module calls from cvmix library use cvmix_tke, only: init_tke, cvmix_coeffs_tke use cvmix_put_get, only: cvmix_put + ! cvmix_kinds_and_types provides cvmix_r8 (CVMix's fixed real64 kind), used below to + ! convert WP<->cvmix_r8 at the CVMix call boundary for single-precision FESOM builds. use cvmix_kinds_and_types use g_cvmix_idemix, only: iwe_n, iwe_Tdis_n, iwe_alpha_c_n @@ -255,21 +257,21 @@ subroutine init_cvmix_tke(partit, mesh) !_______________________________________________________________________ ! call tke initialisation routine from cvmix library - call init_tke(c_k = tke_c_k, & - c_eps = tke_c_eps, & - cd = tke_cd, & - alpha_tke = tke_alpha, & - mxl_min = tke_mxl_min, & - kappaM_min = tke_kappaM_min, & - kappaM_max = tke_kappaM_max, & + call init_tke(c_k = real(tke_c_k, cvmix_r8), & + c_eps = real(tke_c_eps, cvmix_r8), & + cd = real(tke_cd, cvmix_r8), & + alpha_tke = real(tke_alpha, cvmix_r8), & + mxl_min = real(tke_mxl_min, cvmix_r8), & + kappaM_min = real(tke_kappaM_min, cvmix_r8), & + kappaM_max = real(tke_kappaM_max, cvmix_r8), & tke_mxl_choice = tke_mxl_choice, & use_ubound_dirichlet = use_ubound_dirichlet, & use_lbound_dirichlet = use_lbound_dirichlet, & only_tke = tke_only, & l_lc = tke_dolangmuir, & - clc = tke_clangmuir, & - tke_min = tke_min, & - tke_surf_min = tke_surf_min ) + clc = real(tke_clangmuir, cvmix_r8), & + tke_min = real(tke_min, cvmix_r8), & + tke_surf_min = real(tke_surf_min, cvmix_r8) ) end subroutine init_cvmix_tke ! ! @@ -286,6 +288,22 @@ subroutine calc_cvmix_tke(dynamics, partit, mesh) real(kind=WP) :: dz_trr(mesh%nl), bvfreq2(mesh%nl), vshear2(mesh%nl) real(kind=WP) :: tke_Av_old(mesh%nl), tke_Kv_old(mesh%nl), tke_old(mesh%nl) real(kind=WP), dimension(:,:,:), pointer :: UVnode + ! CVMix is a fixed double-precision library (kind cvmix_r8); FESOM's WP may be + ! real4. These buffers convert WP<->cvmix_r8 at the CVMix TKE call boundary so + ! CVMix always runs in double precision. In a DP build cvmix_r8==WP, so the + ! casts/copies are exact no-ops and the result is bit-identical. + ! input arrays (copy-in WP -> cvmix_r8 before the call): + real(cvmix_r8) :: dzw_r8(mesh%nl), dzt_r8(mesh%nl) + real(cvmix_r8) :: tke_old_r8(mesh%nl), av_old_r8(mesh%nl), kv_old_r8(mesh%nl) + real(cvmix_r8) :: ssqr_r8(mesh%nl), nsqr_r8(mesh%nl) + real(cvmix_r8) :: alphac_r8(mesh%nl), eiw_r8(mesh%nl), iwdis_r8(mesh%nl) + real(cvmix_r8) :: plc_r8(mesh%nl) + ! output arrays (copy-out cvmix_r8 -> WP after the call): + real(cvmix_r8) :: tke_new_r8(mesh%nl), kappam_r8(mesh%nl), kappah_r8(mesh%nl) + real(cvmix_r8) :: tbpr_r8(mesh%nl), tspr_r8(mesh%nl), tdif_r8(mesh%nl), tdis_r8(mesh%nl) + real(cvmix_r8) :: twin_r8(mesh%nl), tiwf_r8(mesh%nl), tbck_r8(mesh%nl), ttot_r8(mesh%nl) + real(cvmix_r8) :: lmix_r8(mesh%nl), pr_r8(mesh%nl) + real(cvmix_r8) :: int1_r8(mesh%nl), int2_r8(mesh%nl), int3_r8(mesh%nl) #include "../associate_part_def.h" #include "../associate_mesh_def.h" @@ -433,54 +451,85 @@ subroutine calc_cvmix_tke(dynamics, partit, mesh) tke_Av_old = tke_Av(:,node) tke_Kv_old = tke_Kv(:,node) tke_old = tke(:,node) - + + ! copy-in WP -> cvmix_r8 for the double-precision CVMix TKE call + dzw_r8(nun:nln) = real(hnode(nun:nln,node), cvmix_r8) + dzt_r8(nun:nln+1) = real(dz_trr(nun:nln+1), cvmix_r8) + tke_old_r8(nun:nln+1) = real(tke_old(nun:nln+1), cvmix_r8) + av_old_r8(nun:nln+1) = real(tke_Av_old(nun:nln+1), cvmix_r8) + kv_old_r8(nun:nln+1) = real(tke_Kv_old(nun:nln+1), cvmix_r8) + ssqr_r8(nun:nln+1) = real(vshear2(nun:nln+1), cvmix_r8) + nsqr_r8(nun:nln+1) = real(bvfreq2(nun:nln+1), cvmix_r8) + alphac_r8(nun:nln+1) = real(tke_in3d_iwealphac(nun:nln+1,node), cvmix_r8) + eiw_r8(nun:nln+1) = real(tke_in3d_iwe(nun:nln+1,node), cvmix_r8) + iwdis_r8(nun:nln+1) = real(tke_in3d_iwdis(nun:nln+1,node), cvmix_r8) + plc_r8(nun:nln+1) = real(tke_langmuir(nun:nln+1,node), cvmix_r8) + call cvmix_coeffs_tke(& ! parameter - dzw = hnode(nun:nln,node), & ! distance between layer interface --> hnode - dzt = dz_trr(nun:nln+1), & ! distnace between tracer points + dzw = dzw_r8(nun:nln), & ! distance between layer interface --> hnode + dzt = dzt_r8(nun:nln+1), & ! distnace between tracer points ! nlev = nln, & nlev = nln-nun+1, & max_nlev = nl-1, & - dtime = dt, & - rho_ref = density_0, & - grav = g, & + dtime = real(dt, cvmix_r8), & + rho_ref = real(density_0, cvmix_r8), & + grav = real(g, cvmix_r8), & ! essentials - tke_new = tke( nun:nln+1,node), & ! out--> turbulent kinetic energy - KappaM_out = tke_Av( nun:nln+1,node), & ! out - KappaH_out = tke_Kv( nun:nln+1,node), & ! out - tke_old = tke_old( nun:nln+1), & ! in --> turbulent kinetic energy previous time step - old_KappaM = tke_Av_old(nun:nln+1), & ! in - old_KappaH = tke_Kv_old(nun:nln+1), & ! in - Ssqr = vshear2( nun:nln+1), & ! in --> square vert. vel. shear - Nsqr = bvfreq2( nun:nln+1), & ! in --> square brunt Väisälä freq - alpha_c = tke_in3d_iwealphac(nun:nln+1,node), & ! in for IDEMIX Ri - E_iw = tke_in3d_iwe(nun:nln+1,node), & ! in for IDEMIX Ri + tke_new = tke_new_r8(nun:nln+1), & ! out--> turbulent kinetic energy + KappaM_out = kappam_r8(nun:nln+1), & ! out + KappaH_out = kappah_r8(nun:nln+1), & ! out + tke_old = tke_old_r8(nun:nln+1), & ! in --> turbulent kinetic energy previous time step + old_KappaM = av_old_r8(nun:nln+1), & ! in + old_KappaH = kv_old_r8(nun:nln+1), & ! in + Ssqr = ssqr_r8(nun:nln+1), & ! in --> square vert. vel. shear + Nsqr = nsqr_r8(nun:nln+1), & ! in --> square brunt Väisälä freq + alpha_c = alphac_r8(nun:nln+1), & ! in for IDEMIX Ri + E_iw = eiw_r8(nun:nln+1), & ! in for IDEMIX Ri ! forcing - forc_tke_surf= tke_forc2d_normstress( node), & ! in --> wind stress - forc_rho_surf= tke_forc2d_rhosurf( node), & ! in - bottom_fric = tke_forc2d_botfrict( node), & ! in - iw_diss = tke_in3d_iwdis(nun:nln+1,node), & ! in + forc_tke_surf= real(tke_forc2d_normstress(node), cvmix_r8), & ! in --> wind stress + forc_rho_surf= real(tke_forc2d_rhosurf(node), cvmix_r8), & ! in + bottom_fric = real(tke_forc2d_botfrict(node), cvmix_r8), & ! in + iw_diss = iwdis_r8(nun:nln+1), & ! in ! diagnostics - tke_plc = tke_langmuir(nun:nln+1,node), & ! in - tke_Tbpr = tke_Tbpr(nun:nln+1,node), & ! buoyancy production - tke_Tspr = tke_Tspr(nun:nln+1,node), & ! shear production - tke_Tdif = tke_Tdif(nun:nln+1,node), & ! vertical diffusion d/dz(k d/dz)TKE - tke_Tdis = tke_Tdis(nun:nln+1,node), & ! dissipation - tke_Twin = tke_Twin(nun:nln+1,node), & ! wind forcing - tke_Tiwf = tke_Tiwf(nun:nln+1,node), & ! internal wave forcing when idemix is used - tke_Tbck = tke_Tbck(nun:nln+1,node), & ! background forcing only active if IDEMIX is not active, forcing that results from resetting TKE to minimum background TKE value - tke_Ttot = tke_Ttot(nun:nln+1,node), & ! sum of all terms - tke_Lmix = tke_Lmix(nun:nln+1,node), & ! mixing length scale of the TKE scheme - tke_Pr = tke_Pr( nun:nln+1,node), & ! Prantl number + tke_plc = plc_r8(nun:nln+1), & ! in + tke_Tbpr = tbpr_r8(nun:nln+1), & ! buoyancy production + tke_Tspr = tspr_r8(nun:nln+1), & ! shear production + tke_Tdif = tdif_r8(nun:nln+1), & ! vertical diffusion d/dz(k d/dz)TKE + tke_Tdis = tdis_r8(nun:nln+1), & ! dissipation + tke_Twin = twin_r8(nun:nln+1), & ! wind forcing + tke_Tiwf = tiwf_r8(nun:nln+1), & ! internal wave forcing when idemix is used + tke_Tbck = tbck_r8(nun:nln+1), & ! background forcing only active if IDEMIX is not active, forcing that results from resetting TKE to minimum background TKE value + tke_Ttot = ttot_r8(nun:nln+1), & ! sum of all terms + tke_Lmix = lmix_r8(nun:nln+1), & ! mixing length scale of the TKE scheme + tke_Pr = pr_r8(nun:nln+1), & ! Prantl number ! debugging - cvmix_int_1 = cvmix_dummy_1(nun:nln+1,node), & ! - cvmix_int_2 = cvmix_dummy_2(nun:nln+1,node), & ! - cvmix_int_3 = cvmix_dummy_3(nun:nln+1,node), & ! + cvmix_int_1 = int1_r8(nun:nln+1), & ! + cvmix_int_2 = int2_r8(nun:nln+1), & ! + cvmix_int_3 = int3_r8(nun:nln+1), & ! i = 1, & j = 1, & tstep_count = tstep_count & ) - + + ! copy-out cvmix_r8 -> WP after the double-precision CVMix TKE call + tke(nun:nln+1,node) = real(tke_new_r8(nun:nln+1), WP) + tke_Av(nun:nln+1,node) = real(kappam_r8(nun:nln+1), WP) + tke_Kv(nun:nln+1,node) = real(kappah_r8(nun:nln+1), WP) + tke_Tbpr(nun:nln+1,node) = real(tbpr_r8(nun:nln+1), WP) + tke_Tspr(nun:nln+1,node) = real(tspr_r8(nun:nln+1), WP) + tke_Tdif(nun:nln+1,node) = real(tdif_r8(nun:nln+1), WP) + tke_Tdis(nun:nln+1,node) = real(tdis_r8(nun:nln+1), WP) + tke_Twin(nun:nln+1,node) = real(twin_r8(nun:nln+1), WP) + tke_Tiwf(nun:nln+1,node) = real(tiwf_r8(nun:nln+1), WP) + tke_Tbck(nun:nln+1,node) = real(tbck_r8(nun:nln+1), WP) + tke_Ttot(nun:nln+1,node) = real(ttot_r8(nun:nln+1), WP) + tke_Lmix(nun:nln+1,node) = real(lmix_r8(nun:nln+1), WP) + tke_Pr(nun:nln+1,node) = real(pr_r8(nun:nln+1), WP) + cvmix_dummy_1(nun:nln+1,node) = real(int1_r8(nun:nln+1), WP) + cvmix_dummy_2(nun:nln+1,node) = real(int2_r8(nun:nln+1), WP) + cvmix_dummy_3(nun:nln+1,node) = real(int3_r8(nun:nln+1), WP) + tke_Av(nln+1,node)=0.0_WP tke_Kv(nln+1,node)=0.0_WP tke_Av(nun ,node)=0.0_WP From 76ac7a3b1611b4bdf05ca86e4d8eef5277b7954a Mon Sep 17 00:00:00 2001 From: Suvarchal Kumar Cheedela Date: Sat, 4 Jul 2026 05:48:29 +0200 Subject: [PATCH 5/8] io: report fields requesting output precision above WP (SP awareness) In a single-precision build (WP=4) a namelist.io field may request 8-byte output (e.g. the stock 'salt',1,'m',8). The mean-I/O layer already honours this -- the field is written as NF_DOUBLE with its mean accumulated in real64 -- but the samples are single-precision-sourced. Make the user aware: tally such fields as streams are defined and print one summary line at the end of ini_mean_io (rank 0 only), listing the affected variables. No per-stream or per-step output; silent in double-precision builds where accuracy > WP can never trigger. Add an SP-only ctest (integration_pi_sp_mpi2, label single_precision) that runs the pi mesh and asserts the SINGLE PRECISION MODE startup banner. To support it, add_fesom_test[_with_options] gains optional EXTRA_SUCCESS_MARKERS and LABEL arguments; default behaviour is unchanged. --- cmake/FesomTesting.cmake | 25 +++++++++++++++------ config/namelist.io | 4 ++++ src/io_meandata.F90 | 37 ++++++++++++++++++++++++++++++-- tests/integration/CMakeLists.txt | 15 +++++++++++++ 4 files changed, 73 insertions(+), 8 deletions(-) diff --git a/cmake/FesomTesting.cmake b/cmake/FesomTesting.cmake index e2e72fb83..af6926355 100644 --- a/cmake/FesomTesting.cmake +++ b/cmake/FesomTesting.cmake @@ -366,10 +366,10 @@ endfunction() # Function to add a FESOM integration test with custom options function(add_fesom_test_with_options TEST_NAME MESH_NAME STEP_PER_DAY RUN_LENGTH RUN_LENGTH_UNIT RESTART_LENGTH RESTART_LENGTH_UNIT LOGFILE_OUTFREQ FORCE_ROTATION USE_CAVITY) set(options MPI_TEST) - set(oneValueArgs NP TIMEOUT) - set(multiValueArgs COMMAND_ARGS) + set(oneValueArgs NP TIMEOUT LABEL) + set(multiValueArgs COMMAND_ARGS EXTRA_SUCCESS_MARKERS) cmake_parse_arguments(FESOM_TEST "${options}" "${oneValueArgs}" "${multiValueArgs}" ${ARGN}) - + # Set defaults if(NOT DEFINED FESOM_TEST_NP) set(FESOM_TEST_NP 1) @@ -377,6 +377,14 @@ function(add_fesom_test_with_options TEST_NAME MESH_NAME STEP_PER_DAY RUN_LENGTH if(NOT DEFINED FESOM_TEST_TIMEOUT) set(FESOM_TEST_TIMEOUT 300) # 5 minutes default endif() + + # Assemble the success-marker list: the mandatory clean-exit marker plus any + # caller-supplied EXTRA_SUCCESS_MARKERS (e.g. the single-precision banner). + # Each is quoted so check_fesom_run() receives them as distinct list items. + set(_success_markers "\"fesom should stop with exit status = 0\"") + foreach(_m IN LISTS FESOM_TEST_EXTRA_SUCCESS_MARKERS) + string(APPEND _success_markers " \"${_m}\"") + endforeach() # Create test run directory set(TEST_RUN_DIR "${CMAKE_CURRENT_BINARY_DIR}/${TEST_NAME}") @@ -434,7 +442,7 @@ function(add_fesom_test_with_options TEST_NAME MESH_NAME STEP_PER_DAY RUN_LENGTH RESULT \"\${test_result}\" OUTPUT_LOG \"${TEST_RUN_DIR}/test_output.log\" ERROR_LOG \"${TEST_RUN_DIR}/test_error.log\" - SUCCESS_MARKERS \"fesom should stop with exit status = 0\" + SUCCESS_MARKERS ${_success_markers} REQUIRED_ARTIFACTS \"${RESULT_DIR}/sst.fesom.1948.nc\" ) ") @@ -467,7 +475,7 @@ function(add_fesom_test_with_options TEST_NAME MESH_NAME STEP_PER_DAY RUN_LENGTH RESULT \"\${test_result}\" OUTPUT_LOG \"${TEST_RUN_DIR}/test_output.log\" ERROR_LOG \"${TEST_RUN_DIR}/test_error.log\" - SUCCESS_MARKERS \"fesom should stop with exit status = 0\" + SUCCESS_MARKERS ${_success_markers} REQUIRED_ARTIFACTS \"${RESULT_DIR}/sst.fesom.1948.nc\" ) ") @@ -487,7 +495,12 @@ function(add_fesom_test_with_options TEST_NAME MESH_NAME STEP_PER_DAY RUN_LENGTH TIMEOUT ${FESOM_TEST_TIMEOUT} WORKING_DIRECTORY ${TEST_RUN_DIR} ) - + + # Optional test label (e.g. single_precision) for `ctest -L` + if(DEFINED FESOM_TEST_LABEL) + set_tests_properties(${TEST_NAME} PROPERTIES LABELS "${FESOM_TEST_LABEL}") + endif() + # For MPI tests, set required properties if(FESOM_TEST_MPI_TEST AND FESOM_TEST_NP GREATER 1) set_tests_properties(${TEST_NAME} PROPERTIES diff --git a/config/namelist.io b/config/namelist.io index 27e611ed3..eee53f45d 100644 --- a/config/namelist.io +++ b/config/namelist.io @@ -51,6 +51,10 @@ compression_level = 1 ! compression level for netCDF output (1=fastest, 9 ! frequency = output frequency (integer) ! unit = 'y' (yearly), 'm' (monthly), 'd' (daily), 'h' (hourly), 's' (steps) ! precision = 4 (single precision) or 8 (double precision) +! Note: in a single-precision build (WP=4) a field requesting precision 8 is +! still honoured (written as NF_DOUBLE; means accumulated in real64), but its +! samples are single-precision-sourced. FESOM prints one summary line at +! startup listing such fields. ! ============================================================================ &nml_list io_list = 'sst ',1, 'm', 4, diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index ea805d35a..03e9d8937 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -130,6 +130,13 @@ module io_MEANDATA end type io_entry type(io_entry), save, allocatable, target :: io_list(:) + + ! Fields whose requested output precision exceeds the working precision + ! (e.g. 8-byte output requested in a WP=4 single-precision build). Collected + ! during stream definition and reported once at the end of ini_mean_io so the + ! user is aware, without any per-stream or per-step printing. + integer, save :: n_wp_promoted = 0 + character(len=20), save :: wp_promoted_names(256) ! !-------------------------------------------------------------------------------------------- ! Type for 0D (scalar) output streams - global values with time dimension only @@ -397,7 +404,8 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh) end if end do -!_______________________________________________________________________________ + n_wp_promoted = 0 ! reset the WP-vs-requested-precision awareness tally +!_______________________________________________________________________________ DO i=1, io_listsize SELECT CASE (trim(io_list(i)%id)) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!2D streams!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -2162,7 +2170,22 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh) write(*,*) ' Override via io_list entries in namelist.io' end if end if - + + !___________________________________________________________________________ + ! Awareness: in a single-precision build (WP=4) some output fields may + ! request 8-byte precision in namelist.io. Those are honoured (written as + ! NF_DOUBLE, with means accumulated in real64), but their samples are + ! single-precision-sourced. Report the list once per run (silent in DP, + ! where WP=8 and the condition can never trigger). + if (n_wp_promoted > 0 .and. mype == 0) then + write(*,'(a,i0,a,i0,a)') ' FESOM I/O (WP=', WP, '): ', n_wp_promoted, & + ' field(s) request 8-byte output -> honoured as NF_DOUBLE (real64-accumulated means, SP-sourced samples):' + write(*,'(4x,*(1x,a))') (trim(wp_promoted_names(j)), & + j=1, min(n_wp_promoted, size(wp_promoted_names))) + if (n_wp_promoted > size(wp_promoted_names)) & + write(*,'(4x,a)') '... (list truncated)' + end if + end subroutine @@ -3565,6 +3588,16 @@ subroutine def_stream_after_dimension_specific(entry, name, description, units, !___________________________________________________________________________ entry%accuracy = accuracy + ! Awareness (single-precision builds): note fields whose requested output + ! precision exceeds the working precision (e.g. 8-byte output with WP=4). + ! The request is still honoured (written as NF_DOUBLE, means accumulated in + ! real64), but the samples are single-precision-sourced. Only tallied here; + ! reported once per run at the end of ini_mean_io. + if (accuracy > WP) then + n_wp_promoted = n_wp_promoted + 1 + if (n_wp_promoted <= size(wp_promoted_names)) wp_promoted_names(n_wp_promoted) = name + end if + if (accuracy == i_real8) then allocate(data_strategy_nf_double_type :: entry%data_strategy) elseif (accuracy == i_real4) then diff --git a/tests/integration/CMakeLists.txt b/tests/integration/CMakeLists.txt index a389dfc63..7ae09c80b 100644 --- a/tests/integration/CMakeLists.txt +++ b/tests/integration/CMakeLists.txt @@ -33,6 +33,21 @@ if(ENABLE_MPI_TESTS AND NOT FESOM_COUPLED AND NOT OIFS_COUPLED AND NOT USE_YAC) TIMEOUT ${TEST_TIMEOUT} ) + # Single-precision build: same pi run, but additionally assert that the + # binary is genuinely single precision via the startup banner. This marker + # is only meaningful in an SP build, so the test is only registered there; + # DP builds are unaffected. + if(USE_SINGLE_PRECISION) + add_fesom_test(integration_pi_sp_mpi2 + MPI_TEST + NP 2 + TIMEOUT ${TEST_TIMEOUT} + LABEL single_precision + EXTRA_SUCCESS_MARKERS + "SINGLE PRECISION MODE" + ) + endif() + # Additional MPI test configurations can be added here # For example, testing different domain decompositions elseif(NOT ENABLE_MPI_TESTS) From 5b3d75004975ef7e84317bf958059e58c5c5edc9 Mon Sep 17 00:00:00 2001 From: Suvarchal Kumar Cheedela Date: Sat, 4 Jul 2026 06:11:35 +0200 Subject: [PATCH 6/8] tests: add CVMix (cvmix_TKE, cvmix_KPP) integration tests for DP and SP No integration test exercised the CVMix vertical-mixing path. The stock namelist.oce ships mix_scheme='KPP' -- FESOM's native (non-CVMix) KPP -- so a CVMix-enabled build (CVMIX=ON, the default) never runs CVMix code by default; CI only ever covered native KPP. CVMix schemes such as cvmix_TKE (used at high resolution) went untested. Add integration_pi_cvmix_tke_mpi2 and integration_pi_cvmix_kpp_mpi2 (label cvmix): pi-mesh runs with mix_scheme set to 'cvmix_TKE' / 'cvmix_KPP'. CVMix is a fixed double-precision library (kind cvmix_r8 = real64) built into FESOM's tree via FetchContent; only the FESOM-side WP<->cvmix_r8 conversion differs by build. So each test covers the previously-untested FESOM-DP+CVMix path in a DP build and the FESOM-SP+CVMix path (the conversion shims from the CVMix-in-DP change) in an SP build -- run `ctest -L cvmix` in each build dir. Gated on CVMIX (default ON); no precision guard. Pass/fail is the clean run + artifacts, so a wrong scheme string or a CVMix-less build fails at setup. To support it, add_fesom_test[_with_options] gains an optional MIX_SCHEME argument that rewrites mix_scheme in the test's namelist.oce (whitespace- tolerant, comment-preserving, matching the harness's other rewrites). Verified: both schemes pass in GNU DP and GNU SP builds on pi, with CVMIX_TKE / CVMIX_KPP genuinely initialised (real64 vs real32 params in the logs). --- cmake/FesomTesting.cmake | 15 +++++++++++++-- tests/integration/CMakeLists.txt | 30 ++++++++++++++++++++++++++++++ 2 files changed, 43 insertions(+), 2 deletions(-) diff --git a/cmake/FesomTesting.cmake b/cmake/FesomTesting.cmake index af6926355..6d067c21f 100644 --- a/cmake/FesomTesting.cmake +++ b/cmake/FesomTesting.cmake @@ -366,7 +366,7 @@ endfunction() # Function to add a FESOM integration test with custom options function(add_fesom_test_with_options TEST_NAME MESH_NAME STEP_PER_DAY RUN_LENGTH RUN_LENGTH_UNIT RESTART_LENGTH RESTART_LENGTH_UNIT LOGFILE_OUTFREQ FORCE_ROTATION USE_CAVITY) set(options MPI_TEST) - set(oneValueArgs NP TIMEOUT LABEL) + set(oneValueArgs NP TIMEOUT LABEL MIX_SCHEME) set(multiValueArgs COMMAND_ARGS EXTRA_SUCCESS_MARKERS) cmake_parse_arguments(FESOM_TEST "${options}" "${oneValueArgs}" "${multiValueArgs}" ${ARGN}) @@ -483,7 +483,18 @@ function(add_fesom_test_with_options TEST_NAME MESH_NAME STEP_PER_DAY RUN_LENGTH # Configure namelists for this test with custom options configure_fesom_namelists_with_options("${TEST_RUN_DIR}" "${TEST_DATA_DIR}" "${RESULT_DIR}" "${MESH_NAME}" "${STEP_PER_DAY}" "${RUN_LENGTH}" "${RUN_LENGTH_UNIT}" "${RESTART_LENGTH}" "${RESTART_LENGTH_UNIT}" "${LOGFILE_OUTFREQ}" "${FORCE_ROTATION}" "${USE_CAVITY}") - + + # Optional: override the vertical mixing scheme in namelist.oce (e.g. to + # exercise a CVMix scheme). Whitespace-tolerant and anchored on the char + # before the key, matching the harness's other namelist rewrites; the + # trailing comment is preserved. + if(DEFINED FESOM_TEST_MIX_SCHEME) + file(READ "${TEST_RUN_DIR}/namelist.oce" _oce_content) + string(REGEX REPLACE "([^A-Za-z0-9_]mix_scheme[ \t]*=[ \t]*)'[^']*'" + "\\1'${FESOM_TEST_MIX_SCHEME}'" _oce_content "${_oce_content}") + file(WRITE "${TEST_RUN_DIR}/namelist.oce" "${_oce_content}") + endif() + # Add the test add_test( NAME ${TEST_NAME} diff --git a/tests/integration/CMakeLists.txt b/tests/integration/CMakeLists.txt index 7ae09c80b..9e7660a35 100644 --- a/tests/integration/CMakeLists.txt +++ b/tests/integration/CMakeLists.txt @@ -33,6 +33,36 @@ if(ENABLE_MPI_TESTS AND NOT FESOM_COUPLED AND NOT OIFS_COUPLED AND NOT USE_YAC) TIMEOUT ${TEST_TIMEOUT} ) + # CVMix vertical mixing: run the pi mesh with CVMix schemes so the CVMix code + # path is exercised end to end. NOTE: the stock namelist.oce ships + # mix_scheme='KPP' -- FESOM's native (non-CVMix) KPP -- so a CVMix-enabled + # build (CVMIX=ON) never runs CVMix by default; the scheme must be selected + # explicitly. These tests do that via MIX_SCHEME. + # + # CVMix is a fixed double-precision library (kind cvmix_r8 = real64) built + # into FESOM's tree; the FESOM-side WP<->cvmix_r8 conversion is what differs + # between builds. So each test covers the previously-untested FESOM-DP+CVMix + # path in a DP build and the FESOM-SP+CVMix path (the conversion shims) in an + # SP build. Registered only when CVMix is compiled in (CVMIX=ON, the default). + # Pass/fail is the clean run + artifacts: a wrong scheme string or missing + # CVMix aborts at setup. + if(CVMIX) + add_fesom_test(integration_pi_cvmix_tke_mpi2 + MPI_TEST + NP 2 + TIMEOUT ${TEST_TIMEOUT} + LABEL cvmix + MIX_SCHEME "cvmix_TKE" + ) + add_fesom_test(integration_pi_cvmix_kpp_mpi2 + MPI_TEST + NP 2 + TIMEOUT ${TEST_TIMEOUT} + LABEL cvmix + MIX_SCHEME "cvmix_KPP" + ) + endif() + # Single-precision build: same pi run, but additionally assert that the # binary is genuinely single precision via the startup banner. This marker # is only meaningful in an SP build, so the test is only registered there; From d224cdf95aa46a56821354935cb8a4ce94364f11 Mon Sep 17 00:00:00 2001 From: Suvarchal Kumar Cheedela Date: Mon, 6 Jul 2026 01:37:36 +0200 Subject: [PATCH 7/8] ci/tests: model precision as a build dimension; add SP-only CTest workflow Precision is a property of the build, not the test: in a USE_SINGLE_PRECISION build every integration test runs single precision; in a DP build every test runs double. The previous per-test SP modelling (integration_pi_sp_mpi2 + a single_precision label) was redundant and left the base pi tests unlabelled. Harness (cmake/FesomTesting.cmake): - Every fesom.x integration test now auto-asserts the build's startup banner ('SINGLE'/'DOUBLE PRECISION MODE', chosen by USE_SINGLE_PRECISION) as a required success marker, so each test self-certifies it ran at the intended precision -- a mis-configured build fails loudly on every test. - Every integration test gets a base 'integration' label; the optional LABEL arg (e.g. cvmix) is now appended, so cvmix tests carry 'integration;cvmix'. `ctest -L integration` runs the suite, `ctest -L cvmix` isolates CVMix. tests/integration/CMakeLists.txt: - Remove the now-redundant integration_pi_sp_mpi2 (its guarantee holds for every test in an SP build) and the single_precision label. CI (.github/workflows/fesom2_sp_ctest.yml, new): - SP-only, non-blocking (continue-on-error) workflow: builds full-config in single precision (mirrors DP: MESHPARTITIONER + MESHDIAG ON) so the whole local suite registers in SP, and runs it via `ctest -E '^(remote_|ecbundle_)'`. DP is already covered by fesom2_ctest_runner.yml and the default preset (which already runs the cvmix tests in DP), so no dp/sp matrix -- that would be a third DP run. Promote SP to blocking later by removing continue-on-error; widening to SP meshpart/meshdiag/remote needs no build change, only run selection. Verified: build-dp/build-sp reconfigure; integration_pi_sp_mpi2 gone; 'integration'/'cvmix' labels present, 'single_precision' gone; generated test scripts carry the correct DOUBLE/SINGLE marker; a DP integration test passes with the new required marker. --- .github/workflows/fesom2_sp_ctest.yml | 148 ++++++++++++++++++++++++++ cmake/FesomTesting.cmake | 20 +++- tests/integration/CMakeLists.txt | 19 +--- 3 files changed, 171 insertions(+), 16 deletions(-) create mode 100644 .github/workflows/fesom2_sp_ctest.yml diff --git a/.github/workflows/fesom2_sp_ctest.yml b/.github/workflows/fesom2_sp_ctest.yml new file mode 100644 index 000000000..04c342773 --- /dev/null +++ b/.github/workflows/fesom2_sp_ctest.yml @@ -0,0 +1,148 @@ +name: FESOM2 SP CTests (single precision) + +# Single-precision (USE_SINGLE_PRECISION=ON) counterpart to the DP CTest runner. +# Precision is a property of the BUILD, not the test: this job builds fesom.x (and +# the mesh tools) in single precision with the SAME full config as the DP build, so +# the whole local test suite -- integration (incl. the CVMix schemes) plus the local +# meshpartitioner/meshdiag tests -- registers and runs in SP. Each fesom.x test +# additionally self-asserts the "SINGLE PRECISION MODE" startup banner. +# +# It is OPTIONAL / non-blocking initially (continue-on-error + not a required check): +# DP is already covered by fesom2_ctest_runner.yml and the `default` build preset, so +# this only adds the previously-missing SP coverage. Promote to blocking later by +# removing `continue-on-error`. Widening to SP meshpart/meshdiag/remote needs no build +# change -- only relaxing the ctest selection below. + +on: + pull_request: + types: + - opened + - synchronize + - reopened + branches: [ main ] + + workflow_dispatch: + inputs: + test_timeout: + description: 'Test timeout in seconds' + required: false + default: '600' + type: string + +jobs: + sp_ctests: + runs-on: ubuntu-latest + + # Optional / non-blocking while SP support is being stabilised. Remove this line + # (and add the job to branch-protection required checks) to make SP gate merges. + continue-on-error: true + + env: + CC: gcc + CXX: g++ + FC: gfortran + OMPI_MCA_rmaps_base_oversubscribe: yes + + steps: + - name: Checkout repository + uses: actions/checkout@v4 + + - name: Install system dependencies + run: | + echo "Installing FESOM2 dependencies..." + sudo apt-get update + sudo apt-get install -y \ + gcc \ + gfortran \ + g++ \ + cmake \ + make \ + openmpi-bin \ + libopenmpi-dev \ + libnetcdf-dev \ + libnetcdff-dev \ + pkg-config \ + git \ + wget \ + ca-certificates \ + gnupg \ + lsb-release + + echo "Installing CMake 3.25+ from Kitware repository..." + wget -O - https://apt.kitware.com/keys/kitware-archive-latest.asc 2>/dev/null | gpg --dearmor - | sudo tee /etc/apt/trusted.gpg.d/kitware.gpg >/dev/null + echo "deb https://apt.kitware.com/ubuntu/ $(lsb_release -cs) main" | sudo tee /etc/apt/sources.list.d/kitware.list >/dev/null + sudo apt-get update + sudo apt-get install -y cmake + + echo "Verifying installations..." + gcc --version + gfortran --version + mpirun --version + cmake --version + pkg-config --modversion netcdf + pkg-config --modversion netcdf-fortran + + - name: Set up environment variables + run: | + echo "NETCDF_ROOT=/usr" >> $GITHUB_ENV + echo "NETCDF_Fortran_ROOT=/usr" >> $GITHUB_ENV + echo "FCFLAGS=-fallow-argument-mismatch" >> $GITHUB_ENV + echo "FFLAGS=-fallow-argument-mismatch" >> $GITHUB_ENV + echo "OMPI_MCA_btl_vader_single_copy_mechanism=none" >> $GITHUB_ENV + echo "OMPI_MCA_rmaps_base_oversubscribe=yes" >> $GITHUB_ENV + echo "OMPI_MCA_btl_base_warn_component_unused=0" >> $GITHUB_ENV + + - name: Configure FESOM2 (single precision, full config = mirrors DP) + run: | + echo "Configuring FESOM2 single-precision build..." + mkdir -p build + cd build + cmake \ + -DCMAKE_BUILD_TYPE=Release \ + -DUSE_SINGLE_PRECISION=ON \ + -DBUILD_TESTING=ON \ + -DENABLE_MPI_TESTS=ON \ + -DBUILD_MESHDIAG=ON \ + -DBUILD_MESHPARTITIONER=ON \ + -DTEST_TIMEOUT=${{ github.event.inputs.test_timeout || '600' }} \ + -DCMAKE_C_COMPILER=gcc \ + -DCMAKE_CXX_COMPILER=g++ \ + -DCMAKE_Fortran_COMPILER=gfortran \ + -DMPI_C_COMPILER=mpicc \ + -DMPI_CXX_COMPILER=mpicxx \ + -DMPI_Fortran_COMPILER=mpifort \ + -DNetCDF_ROOT=/usr \ + -DNetCDF_Fortran_ROOT=/usr \ + -DCMAKE_Fortran_FLAGS="-fallow-argument-mismatch" \ + .. + + - name: Build FESOM2 + run: | + cd build + make -j$(nproc) VERBOSE=1 + if [ -f "bin/fesom.x" ]; then + echo "✅ single-precision fesom.x created" + ls -la bin/fesom.x + else + echo "❌ fesom.x not found" + exit 1 + fi + + - name: List available tests + run: | + cd build + echo "=== Available Tests ===" + ctest -N + echo "=== Labels ===" + ctest --print-labels + + - name: Run CTest (local suite, single precision) + run: | + cd build + echo "Running local tests in single precision (excluding remote + ecbundle)..." + # Local set only: integration_pi_* (incl. cvmix) + local meshpartitioner/meshdiag. + # Remote (network downloads) and ecbundle are excluded. The -E exclusion auto- + # includes any new local test without editing this line. + ctest -E '^(remote_|ecbundle_)' \ + --output-on-failure \ + --timeout ${{ github.event.inputs.test_timeout || '600' }} diff --git a/cmake/FesomTesting.cmake b/cmake/FesomTesting.cmake index 6d067c21f..15817e122 100644 --- a/cmake/FesomTesting.cmake +++ b/cmake/FesomTesting.cmake @@ -385,6 +385,17 @@ function(add_fesom_test_with_options TEST_NAME MESH_NAME STEP_PER_DAY RUN_LENGTH foreach(_m IN LISTS FESOM_TEST_EXTRA_SUCCESS_MARKERS) string(APPEND _success_markers " \"${_m}\"") endforeach() + + # Self-certify the build precision: every fesom.x run prints a working-precision + # banner at startup. Appending the build's banner as a mandatory marker makes each + # integration test verify it actually ran at the intended precision -- so a + # mis-configured build (DP when SP was wanted, or vice versa) fails loudly on every + # test rather than silently running the wrong precision. Silent no-op in DP. + if(USE_SINGLE_PRECISION) + string(APPEND _success_markers " \"SINGLE PRECISION MODE\"") + else() + string(APPEND _success_markers " \"DOUBLE PRECISION MODE\"") + endif() # Create test run directory set(TEST_RUN_DIR "${CMAKE_CURRENT_BINARY_DIR}/${TEST_NAME}") @@ -507,10 +518,15 @@ function(add_fesom_test_with_options TEST_NAME MESH_NAME STEP_PER_DAY RUN_LENGTH WORKING_DIRECTORY ${TEST_RUN_DIR} ) - # Optional test label (e.g. single_precision) for `ctest -L` + # Every fesom.x integration test carries the base label 'integration', plus any + # extra label passed via LABEL (e.g. cvmix). CTest LABELS is a list, so this yields + # e.g. "integration" for the base pi tests and "integration;cvmix" for CVMix tests: + # `ctest -L integration` runs the whole suite, `ctest -L cvmix` isolates CVMix. + set(_labels "integration") if(DEFINED FESOM_TEST_LABEL) - set_tests_properties(${TEST_NAME} PROPERTIES LABELS "${FESOM_TEST_LABEL}") + list(APPEND _labels "${FESOM_TEST_LABEL}") endif() + set_tests_properties(${TEST_NAME} PROPERTIES LABELS "${_labels}") # For MPI tests, set required properties if(FESOM_TEST_MPI_TEST AND FESOM_TEST_NP GREATER 1) diff --git a/tests/integration/CMakeLists.txt b/tests/integration/CMakeLists.txt index 9e7660a35..301eb7fba 100644 --- a/tests/integration/CMakeLists.txt +++ b/tests/integration/CMakeLists.txt @@ -63,20 +63,11 @@ if(ENABLE_MPI_TESTS AND NOT FESOM_COUPLED AND NOT OIFS_COUPLED AND NOT USE_YAC) ) endif() - # Single-precision build: same pi run, but additionally assert that the - # binary is genuinely single precision via the startup banner. This marker - # is only meaningful in an SP build, so the test is only registered there; - # DP builds are unaffected. - if(USE_SINGLE_PRECISION) - add_fesom_test(integration_pi_sp_mpi2 - MPI_TEST - NP 2 - TIMEOUT ${TEST_TIMEOUT} - LABEL single_precision - EXTRA_SUCCESS_MARKERS - "SINGLE PRECISION MODE" - ) - endif() + # Note: there is no separate "SP" integration test. Precision is a property of + # the build, not the test -- in a USE_SINGLE_PRECISION build every test above runs + # single precision, and each fesom.x test auto-asserts the build's precision banner + # (see add_fesom_test_with_options in cmake/FesomTesting.cmake). Run the suite in an + # SP build dir to get SP coverage of everything, including the CVMix schemes. # Additional MPI test configurations can be added here # For example, testing different domain decompositions From 636a60a3523b82e22c83bc679948c67d8e61c372 Mon Sep 17 00:00:00 2001 From: Suvarchal Kumar Cheedela Date: Mon, 6 Jul 2026 02:33:43 +0200 Subject: [PATCH 8/8] ci: fix invalid YAML in fesom2_intel_tests.yml (stray leading dash) Line 51 began with a literal '-' (a diff fragment committed by mistake in #938), which terminated the `run: |` block scalar and made the whole workflow invalid YAML -- GitHub rejects it ("Invalid workflow file") and act cannot parse it. Restore the intended 8-space indentation so the line stays a comment inside the run block; the commented-out ctest invocation is preserved, not deleted. --- .github/workflows/fesom2_intel_tests.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/fesom2_intel_tests.yml b/.github/workflows/fesom2_intel_tests.yml index 812a38692..2739165ad 100644 --- a/.github/workflows/fesom2_intel_tests.yml +++ b/.github/workflows/fesom2_intel_tests.yml @@ -48,5 +48,5 @@ jobs: echo "Intel build OK (integration tests skipped)" # TODO issues with netcdf: cd build -- # ctest -R "integration_*" --output-on-failure --verbose + # ctest -R "integration_*" --output-on-failure --verbose