From 4c3e33a474eabff87e94f367cde111fea1a2460f Mon Sep 17 00:00:00 2001 From: Matt Shin Date: Wed, 27 Jul 2022 14:10:15 +0100 Subject: [PATCH] Import @6240 um13.0 --- Makefile | 2 +- common/src/shumlib_version.h | 2 +- scripts/meto_install_shumlib.sh | 34 +- shum_byteswap/src/Makefile | 8 +- shum_byteswap/src/c_shum_byteswap.c | 22 + shum_byteswap/src/c_shum_byteswap_opt.h | 64 +- .../src/f_shum_horizontal_field_interp.f90 | 601 ++++++++++-------- .../src/f_shum_spiral_search.f90 | 162 +++-- 8 files changed, 508 insertions(+), 387 deletions(-) diff --git a/Makefile b/Makefile index e9c7a34..c9674b6 100644 --- a/Makefile +++ b/Makefile @@ -156,7 +156,7 @@ PACK_PREREQ=STR_CONV CONSTS # Horizontal Interpolation #-------------- HORIZ_INTERP=shum_horizontal_field_interp -HORIZ_INTERP_PREREQ= +HORIZ_INTERP_PREREQ=NUMTOOLS # Thread Utils #-------------- diff --git a/common/src/shumlib_version.h b/common/src/shumlib_version.h index 7029256..b67a32b 100644 --- a/common/src/shumlib_version.h +++ b/common/src/shumlib_version.h @@ -41,7 +41,7 @@ * where "X" is the release number in month "MM" of year "YYYY" */ #if !defined(SHUMLIB_VERSION) -#define SHUMLIB_VERSION 2021101 +#define SHUMLIB_VERSION 2022071 #endif /* 2-stage expansion which will replace in the including code: diff --git a/scripts/meto_install_shumlib.sh b/scripts/meto_install_shumlib.sh index b443f33..060e6fa 100755 --- a/scripts/meto_install_shumlib.sh +++ b/scripts/meto_install_shumlib.sh @@ -27,8 +27,8 @@ # USAGE: (Note - must be run from the toplevel Shumlib directory!) # scripts/meto_install_shumlib.sh [xc40|x86|ex1a] # -# This script was used to install shumlib version 2022.02.1 -# and was intended for use with the UM at UM 12.2 +# This script was used to install shumlib version 2022.07.1 +# and was intended for use with the UM at UM 13.0 # set -eu @@ -644,23 +644,15 @@ if [ $PLATFORM == "xc40" ] || [ $PLATFORM == $THIS ] ; then fi fi -# Crayftn/CrayCC 12.0.3 -THIS="ex1a_cray_12.0.3" +# Crayftn/CrayCC 14.0.0 +THIS="ex1a_cray_14.0.0" if [ $PLATFORM == "ex1a" ] || [ $PLATFORM == $THIS ] ; then ( - module purge - module load PrgEnv-cray/8.2.0 - module switch cce cce/12.0.3 - module switch cray-mpich cray-mpich/8.1.10 - module switch craype craype/2.7.11 - module load craype-x86-milan - module load perftools-base - module load xpmem - module load craype-network-ofi - module load craype-hugepages4M + module switch PrgEnv-cray PrgEnv-cray/8.3.3 + module load cpe/22.05 CONFIG=meto-ex1a-crayftn12.0.1+-craycc - LIBDIR=$BUILD_DESTINATION/meto-ex1a-crayftn-12.0.3-craycc-12.0.3 + LIBDIR=$BUILD_DESTINATION/meto-ex1a-crayftn-14.0.0-craycc-14.0.0 build_openmp_onoff $CONFIG $LIBDIR $(sed -e "s/\bshum_fieldsfile_class\b//g" \ -e "s/\bshum_fieldsfile\b//g" \ <<< $LIB_DIRS) @@ -676,15 +668,9 @@ THIS="ex1a_gnu_10.3.0" if [ $PLATFORM == "ex1a" ] || [ $PLATFORM == $THIS ] ; then ( - module purge - module load PrgEnv-gnu/8.2.0 - module switch gcc gcc/10.3.0 - module load craype-x86-milan - module switch cray-mpich/8.1.10 - module switch craype/2.7.11 - module load perftools-base - module load xpmem - module load craype-network-ofi + module switch PrgEnv-cray PrgEnv-gnu/8.3.3 + module load cpe/22.05 + module switch gcc gcc/10.3.0 CONFIG=meto-ex1a-gfortran-gcc LIBDIR=$BUILD_DESTINATION/meto-ex1a-gfortran-10.3.0-gcc-10.3.0 build_openmp_onoff $CONFIG $LIBDIR all_libs diff --git a/shum_byteswap/src/Makefile b/shum_byteswap/src/Makefile index ae5957f..925b229 100644 --- a/shum_byteswap/src/Makefile +++ b/shum_byteswap/src/Makefile @@ -17,8 +17,8 @@ include ${COMMON_DIR}/Makefile-version # Static library - uses archiver to create library #------------------------------------------------------------------------------- -c_shum_byteswap.o: c_shum_byteswap.c c_shum_byteswap.h c_shum_byteswap_opt.h - ${CC} ${CCFLAGS} -c -I${LIBDIR_OUT}/include $< +c_shum_byteswap.o: c_shum_byteswap.c c_shum_byteswap.h c_shum_byteswap_opt.h ${COMMON_DIR}/c_shum_compile_diag_suspend.h + ${CC} ${CCFLAGS} -c -I${LIBDIR_OUT}/include -I${COMMON_DIR} $< f_shum_byteswap.o: f_shum_byteswap.f90 ${FC} ${FCFLAGS} -c -I${LIBDIR_OUT}/include $< @@ -30,8 +30,8 @@ libshum_byteswap.a: c_shum_byteswap.o f_shum_byteswap.o ${VERSION_OBJECTS} # are suffixed to keep them separate to the above (since the dynamic library # required position-independent-code flags whilst the static version does not) #------------------------------------------------------------------------------- -c_shum_byteswap_PIC.o: c_shum_byteswap.c c_shum_byteswap.h c_shum_byteswap_opt.h - ${CC} ${CCFLAGS} ${CCFLAGS_PIC} -c -I${LIBDIR_OUT}/include $< -o $@ +c_shum_byteswap_PIC.o: c_shum_byteswap.c c_shum_byteswap.h c_shum_byteswap_opt.h ${COMMON_DIR}/c_shum_compile_diag_suspend.h + ${CC} ${CCFLAGS} ${CCFLAGS_PIC} -c -I${LIBDIR_OUT}/include -I${COMMON_DIR} $< -o $@ f_shum_byteswap_PIC.o: f_shum_byteswap.f90 ${FC} ${FCFLAGS} ${FCFLAGS_PIC} -c -I${LIBDIR_OUT}/include $< -o $@ diff --git a/shum_byteswap/src/c_shum_byteswap.c b/shum_byteswap/src/c_shum_byteswap.c index cc1af89..14da368 100644 --- a/shum_byteswap/src/c_shum_byteswap.c +++ b/shum_byteswap/src/c_shum_byteswap.c @@ -31,6 +31,7 @@ #include #include "c_shum_byteswap.h" #include "c_shum_byteswap_opt.h" +#include "c_shum_compile_diag_suspend.h" #if defined(_OPENMP) && defined(SHUM_USE_C_OPENMP_VIA_THREAD_UTILS) #include "c_shum_thread_utils.h" @@ -53,6 +54,21 @@ #define INLINEQUAL static inline #endif +/* We need to suspend compiler diagnostics for the expansion of byteswap macros + * on some systems, as they bring in reserved identifiers from system headers. + */ +#if defined(__clang__) +#if __clang_major__ >= 13 +#define SUSPEND_RES_IDNET SHUM_COMPILE_DIAG_SUSPEND(-Wreserved-identifier) +#define RESUME_RES_IDENT SHUM_COMPILE_DIAG_RESUME +#endif +#endif + +#if !defined(SUSPEND_RES_IDNET) +#define SUSPEND_RES_IDNET +#define RESUME_RES_IDENT +#endif + /* -------------------------------------------------------------------------- */ /* Prototypes and Typedefs */ /* -------------------------------------------------------------------------- */ @@ -169,7 +185,9 @@ INLINEQUAL void c_shum_byteswap_par_swap64(void **const array, #endif for (i=0; i<=span; i=i+1) { + SUSPEND_RES_IDNET ptr_64[min+i*cincr] = (uint64_t)bswap_64(ptr_64[min+i*cincr]); + RESUME_RES_IDENT } } } @@ -193,7 +211,9 @@ INLINEQUAL void c_shum_byteswap_par_swap32(void **const array, #endif for (i=0; i<=span; i=i+1) { + SUSPEND_RES_IDNET ptr_32[min+i*cincr] = (uint32_t)bswap_32(ptr_32[min+i*cincr]); + RESUME_RES_IDENT } } } @@ -217,7 +237,9 @@ INLINEQUAL void c_shum_byteswap_par_swap16(void **const array, #endif for (i=0; i<=span; i=i+1) { + SUSPEND_RES_IDNET ptr_16[min+i*cincr] = (uint16_t)bswap_16(ptr_16[min+i*cincr]); + RESUME_RES_IDENT } } } diff --git a/shum_byteswap/src/c_shum_byteswap_opt.h b/shum_byteswap/src/c_shum_byteswap_opt.h index 6416b39..0db2dee 100644 --- a/shum_byteswap/src/c_shum_byteswap_opt.h +++ b/shum_byteswap/src/c_shum_byteswap_opt.h @@ -64,16 +64,23 @@ * standard C99, so will always work. * * vi) If the user defines one of the macros C_USE_BSWAP_BITSHIFT, - * C_USE_BSWAP_BUILTINS, or C_USE_BSWAP_BYTESWAP_H, override the above - * cases and use the direct bit-shift based macros, the __builtin_bswap*() - * non-standard intrinsic functions, or the headers - * repectively, as per that choice. + * C_USE_BSWAP_BUILTINS, C_USE_BSWAP_BYTESWAP_H, or + * C_USE_BSWAP_OSBYTEORDER_H override the above cases and use the direct + * bit-shift based macros, the __builtin_bswap*() non-standard intrinsic + * functions, the headers, or the + * headers repectively, as per that choice. + * + * Note: the following cases were added to support more systems at a later date + * and so are out of sequence. A user choice (case vi) will still override + * these defaults * * vii) If we are usining a version of the Cray compiler (defines _CRAYC) with * GNU extensions enabled (-hgnu) and has __GNUC__ greater than version * 6.1, fall-back to the non-standard intrinsic __builtin_bswap*() * functions for all data sizes. * + * viii) If we are on a Mac OS X / Darwin system (defines __APPLE__) instead + * use the non-standard header */ /*----------------*/ @@ -87,21 +94,21 @@ #define C_SHUM_BSWAP_HASGNU 0 #endif -#if (C_SHUM_BSWAP_HASGNU && ((__GNUC__ > 4) || \ +#if (C_SHUM_BSWAP_HASGNU && ((__GNUC__ > 4) || \ (__GNUC__ == 4 && __GNUC_MINOR__ >= 3))) #define C_SHUM_BSWAP_HASGNU4_3 1 #else #define C_SHUM_BSWAP_HASGNU4_3 0 #endif -#if (C_SHUM_BSWAP_HASGNU && ((__GNUC__ > 4) || \ +#if (C_SHUM_BSWAP_HASGNU && ((__GNUC__ > 4) || \ (__GNUC__ == 4 && __GNUC_MINOR__ >= 8))) #define C_SHUM_BSWAP_HASGNU4_8 1 #else #define C_SHUM_BSWAP_HASGNU4_8 0 #endif -#if (C_SHUM_BSWAP_HASGNU && ((__GNUC__ > 6) || \ +#if (C_SHUM_BSWAP_HASGNU && ((__GNUC__ > 6) || \ (__GNUC__ == 6 && __GNUC_MINOR__ >= 1))) #define C_SHUM_BSWAP_HASGNU6_1 1 #else @@ -123,27 +130,51 @@ /* test for user choice */ -#if defined C_USE_BSWAP_BYTESWAP_H +#if defined(C_USE_BSWAP_BYTESWAP_H) /* user choice is */ #define C_USE_BSWAP_USERCHOICE +#undef C_USE_BSWAP_OSBYTEORDER_H +#undef C_USE_BSWAP_BITSHIFT +#undef C_USE_BSWAP_BUILTINS + +#elif defined(C_USE_BSWAP_OSBYTEORDER_H) + +/* user choice is */ +#define C_USE_BSWAP_USERCHOICE +#undef C_USE_BSWAP_BYTESWAP_H +#undef C_USE_BSWAP_BITSHIFT +#undef C_USE_BSWAP_BUILTINS -#elif defined C_USE_BSWAP_BUILTINS +#elif defined(C_USE_BSWAP_BUILTINS) /* user choice is __builtin_bswap*() */ #define C_USE_BSWAP_USERCHOICE #define C_USE_BSWAP_BUILTINS_64 #define C_USE_BSWAP_BUILTINS_32 #define C_USE_BSWAP_BUILTINS_16 +#undef C_USE_BSWAP_OSBYTEORDER_H +#undef C_USE_BSWAP_BYTESWAP_H +#undef C_USE_BSWAP_BITSHIFT -#elif defined C_USE_BSWAP_BITSHIFT +#elif defined(C_USE_BSWAP_BITSHIFT) /* user choice is bit-shift macros */ #define C_USE_BSWAP_USERCHOICE #define C_USE_BSWAP_64_BITSHIFT #define C_USE_BSWAP_32_BITSHIFT #define C_USE_BSWAP_16_BITSHIFT +#undef C_USE_BSWAP_OSBYTEORDER_H +#undef C_USE_BSWAP_BYTESWAP_H +#undef C_USE_BSWAP_BUILTINS + +#endif +#if !defined(C_USE_BSWAP_USERCHOICE) +#undef C_USE_BSWAP_OSBYTEORDER_H +#undef C_USE_BSWAP_BYTESWAP_H +#undef C_USE_BSWAP_BITSHIFT +#undef C_USE_BSWAP_BUILTINS #endif /* test for used case */ @@ -152,6 +183,11 @@ /* case vi) */ +#elif defined(__APPLE__) + +/* case viii) */ +#define C_USE_BSWAP_OSBYTEORDER_H + #elif C_SHUM_BSWAP_HASGNU && !C_SHUM_BSWAP_HASGNU4_3 /* case i) */ @@ -214,6 +250,14 @@ #include #endif +#if defined(C_USE_BSWAP_OSBYTEORDER_H) +#include +#define bswap_16(x) OSSwapInt16(x) +#define bswap_32(x) OSSwapInt32(x) +#define bswap_64(x) OSSwapInt64(x) +#endif + + #if defined(C_USE_BSWAP_64_BITSHIFT) #define bswap_64(x) \ diff --git a/shum_horizontal_field_interp/src/f_shum_horizontal_field_interp.f90 b/shum_horizontal_field_interp/src/f_shum_horizontal_field_interp.f90 index f87a2cd..6c5597c 100644 --- a/shum_horizontal_field_interp/src/f_shum_horizontal_field_interp.f90 +++ b/shum_horizontal_field_interp/src/f_shum_horizontal_field_interp.f90 @@ -34,12 +34,12 @@ MODULE f_shum_horizontal_field_interp_mod PRIVATE -PUBLIC :: f_shum_horizontal_field_bi_lin_interp_get_coeffs & - , f_shum_horizontal_field_bi_lin_interp_calc & - , f_shum_find_source_box_indices & - , f_shum_calc_weights & - , f_shum_cart_horizontal_field_bi_lin_interp_get_coeffs & - , f_shum_find_source_cart_box_indices & +PUBLIC :: f_shum_horizontal_field_bi_lin_interp_get_coeffs & + , f_shum_horizontal_field_bi_lin_interp_calc & + , f_shum_find_source_box_indices & + , f_shum_calc_weights & + , f_shum_cart_horizontal_field_bi_lin_interp_get_coeffs & + , f_shum_find_source_cart_box_indices & , f_shum_calc_cart_weights @@ -51,11 +51,17 @@ MODULE f_shum_horizontal_field_interp_mod ! Additional protection for the case that FLOAT/DOUBLE do not conform to the ! ! sizes we expect is provided via the "precision_bomb" macro-file ! !------------------------------------------------------------------------------! - INTEGER, PARAMETER :: int64 = C_INT64_T - INTEGER, PARAMETER :: int32 = C_INT32_T - INTEGER, PARAMETER :: real64 = C_DOUBLE - INTEGER, PARAMETER :: real32 = C_FLOAT - INTEGER, PARAMETER :: bool = C_BOOL +INTEGER, PARAMETER :: shum_int64 = C_INT64_T +INTEGER, PARAMETER :: shum_int32 = C_INT32_T +INTEGER, PARAMETER :: shum_real64 = C_DOUBLE +INTEGER, PARAMETER :: shum_real32 = C_FLOAT +INTEGER, PARAMETER :: shum_bool = C_BOOL +!------------------------------------------------------------------------------! + +INTEGER(KIND=shum_int64), PARAMETER :: zero64 = 0_shum_int64 +INTEGER(KIND=shum_int64), PARAMETER :: one64 = 1_shum_int64 +INTEGER(KIND=shum_int64), PARAMETER :: two64 = 2_shum_int64 + !------------------------------------------------------------------------------! CONTAINS @@ -68,6 +74,8 @@ SUBROUTINE f_shum_horizontal_field_bi_lin_interp_calc & , weight_b_l, weight_b_r, weight_t_l, weight_t_r & , data_out ) +USE f_shum_is_denormal_mod, ONLY: f_shum_is_denormal + IMPLICIT NONE ! Description: @@ -77,62 +85,130 @@ SUBROUTINE f_shum_horizontal_field_bi_lin_interp_calc & ! Subroutine arguments ! Scalar arguments -INTEGER(KIND=int64), INTENT(IN) :: rows_in +INTEGER(KIND=shum_int64), INTENT(IN) :: rows_in ! Number of P rows on source grid -INTEGER(KIND=int64), INTENT(IN) :: row_length_in +INTEGER(KIND=shum_int64), INTENT(IN) :: row_length_in ! Number of pts per row on source grid -INTEGER(KIND=int64), INTENT(IN) :: len_field_out +INTEGER(KIND=shum_int64), INTENT(IN) :: len_field_out ! Number of points on target grid ! Array arguments -INTEGER(KIND=int64), INTENT(IN) :: index_b_l(len_field_out) +INTEGER(KIND=shum_int64), INTENT(IN) :: index_b_l(len_field_out) ! Index of bottom left ! corner of source gridbox -INTEGER(KIND=int64), INTENT(IN) :: index_b_r(len_field_out) +INTEGER(KIND=shum_int64), INTENT(IN) :: index_b_r(len_field_out) ! Index of bottom right ! corner of source gridbox -INTEGER(KIND=int64), INTENT(IN) :: index_t_l(len_field_out) +INTEGER(KIND=shum_int64), INTENT(IN) :: index_t_l(len_field_out) ! Index of top left ! corner of source gridbox -INTEGER(KIND=int64), INTENT(IN) :: index_t_r(len_field_out) +INTEGER(KIND=shum_int64), INTENT(IN) :: index_t_r(len_field_out) ! Index of top right ! corner of source gridbox -REAL(KIND=real64), INTENT(IN) :: data_in(rows_in*row_length_in) +REAL(KIND=shum_real64), INTENT(IN) :: data_in(rows_in*row_length_in) ! Data before interpolation -REAL(KIND=real64), INTENT(IN) :: weight_b_l(len_field_out) +REAL(KIND=shum_real64), INTENT(IN) :: weight_b_l(len_field_out) ! Weight applied to value at bot. ! left corner of source gridbox -REAL(KIND=real64), INTENT(IN) :: weight_b_r(len_field_out) +REAL(KIND=shum_real64), INTENT(IN) :: weight_b_r(len_field_out) ! Weight applied to value at bot. ! right corner of source gridbox -REAL(KIND=real64), INTENT(IN) :: weight_t_l(len_field_out) +REAL(KIND=shum_real64), INTENT(IN) :: weight_t_l(len_field_out) ! Weight applied to value at top ! left corner of source gridbox -REAL(KIND=real64), INTENT(IN) :: weight_t_r(len_field_out) +REAL(KIND=shum_real64), INTENT(IN) :: weight_t_r(len_field_out) ! Weight applied to value at top ! right corner of source gridbox -REAL(KIND=real64), INTENT(OUT) :: data_out(len_field_out) +REAL(KIND=shum_real64), INTENT(OUT) :: data_out(len_field_out) ! Data after interpolation ! Local scalars: -INTEGER(KIND=int64) :: i ! loop index +INTEGER(KIND=shum_int64) :: i ! loop index + +! Local reals +REAL(KIND=shum_real64) :: & + data_out_a, & + data_out_b, & + data_out_c, & + data_out_d + +! Data after interpolation ! 1. Carry out horizontal interpolation !$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) & -!$OMP& SHARED(len_field_out, data_out, weight_b_l, weight_b_r, & -!$OMP& weight_t_l, weight_t_r, index_b_l, index_b_r, & -!$OMP& index_t_l, index_t_r, & -!$OMP& data_in) & -!$OMP& PRIVATE(i) +!$OMP SHARED(len_field_out, data_out, weight_b_l, weight_b_r, & +!$OMP weight_t_l, weight_t_r, index_b_l, index_b_r, & +!$OMP index_t_l, index_t_r, & +!$OMP data_in) & +!$OMP PRIVATE(i, data_out_a, data_out_b, data_out_c, data_out_d ) DO i=1, len_field_out - data_out(i) = weight_b_l(i)*data_in(index_b_l(i)) & - + weight_b_r(i)*data_in(index_b_r(i)) & - + weight_t_l(i)*data_in(index_t_l(i)) & - + weight_t_r(i)*data_in(index_t_r(i)) + ! Weights should be in the range 0.0 <= weight <= 1.0 + ! This means the result cannot overflow. We still need to check for + ! underflow. + + ! b_l component + + IF (f_shum_is_denormal(data_in(index_b_l(i))) .OR. & + f_shum_is_denormal(weight_b_l(i)) .OR. & + weight_b_l(i) < TINY(data_out_a) ) THEN + data_out_a = 0.0 + ELSE + IF (ABS(data_in(index_b_l(i))) < TINY(data_out_a) / weight_b_l(i)) THEN + data_out_a = 0.0 + ELSE + data_out_a = weight_b_l(i)*data_in(index_b_l(i)) + END IF + END IF + + ! b_r component + + IF (f_shum_is_denormal(data_in(index_b_r(i))) .OR. & + f_shum_is_denormal(weight_b_r(i)) .OR. & + weight_b_r(i) < TINY(data_out_b) ) THEN + data_out_b = 0.0 + ELSE + IF (ABS(data_in(index_b_r(i))) < TINY(data_out_b) / weight_b_r(i)) THEN + data_out_b = 0.0 + ELSE + data_out_b = weight_b_r(i)*data_in(index_b_r(i)) + END IF + END IF + + ! t_l component + + IF (f_shum_is_denormal(data_in(index_t_l(i))) .OR. & + f_shum_is_denormal(weight_t_l(i)) .OR. & + weight_t_l(i) < TINY(data_out_c) ) THEN + data_out_c = 0.0 + ELSE + IF (ABS(data_in(index_t_l(i))) < TINY(data_out_c) / weight_t_l(i)) THEN + data_out_c = 0.0 + ELSE + data_out_c = weight_t_l(i)*data_in(index_t_l(i)) + END IF + END IF + + ! t_r component + + IF (f_shum_is_denormal(data_in(index_t_r(i))) .OR. & + f_shum_is_denormal(weight_t_r(i)) .OR. & + weight_t_r(i) < TINY(data_out_d) ) THEN + data_out_d = 0.0 + ELSE + IF (ABS(data_in(index_t_r(i))) < TINY(data_out_d) / weight_t_r(i)) THEN + data_out_d = 0.0 + ELSE + data_out_d = weight_t_r(i)*data_in(index_t_r(i)) + END IF + END IF + + ! combined components + + data_out(i) = data_out_a + data_out_b + data_out_c + data_out_d END DO !$OMP END PARALLEL DO @@ -176,58 +252,59 @@ SUBROUTINE f_shum_horizontal_field_bi_lin_interp_get_coeffs & IMPLICIT NONE -INTEGER(KIND=int64), INTENT(IN) :: points_lambda_srce +INTEGER(KIND=shum_int64), INTENT(IN) :: points_lambda_srce ! Number of lambda points on source grid -INTEGER(KIND=int64), INTENT(IN) :: points_phi_srce +INTEGER(KIND=shum_int64), INTENT(IN) :: points_phi_srce ! Number of phi points on source grid -INTEGER(KIND=int64), INTENT(IN) :: points +INTEGER(KIND=shum_int64), INTENT(IN) :: points ! Total number of points on target grid -INTEGER(KIND=int64), INTENT(OUT) :: index_b_l(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: index_b_l(points) ! Index of bottom left corner ! of source gridbox -INTEGER(KIND=int64), INTENT(OUT) :: index_b_r(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: index_b_r(points) ! Index of bottom right corner ! of source gridbox -INTEGER(KIND=int64), INTENT(OUT) :: index_t_l(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: index_t_l(points) ! Index of top left corner ! of source gridbox -INTEGER(KIND=int64), INTENT(OUT) :: index_t_r(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: index_t_r(points) ! Index of top right corner ! of source gridbox -REAL(KIND=real64), INTENT(IN) :: lambda_targ(points) +REAL(KIND=shum_real64), INTENT(IN) :: lambda_targ(points) ! Lambda coords of target grid in degrees ! using same rotation as source grid -REAL(KIND=real64), INTENT(IN) :: phi_targ(points) +REAL(KIND=shum_real64), INTENT(IN) :: phi_targ(points) ! Phi coords of target grid in degrees using ! same rotation as source grid -REAL(KIND=real64), INTENT(IN) :: lambda_srce(points_lambda_srce) +REAL(KIND=shum_real64), INTENT(IN) :: lambda_srce(points_lambda_srce) ! Lambda coords of source grid in degrees -REAL(KIND=real64), INTENT(IN) :: phi_srce(points_phi_srce) +REAL(KIND=shum_real64), INTENT(IN) :: phi_srce(points_phi_srce) ! Phi coords of source grid in degrees -REAL(KIND=real64), INTENT(OUT) :: weight_t_r(points) +REAL(KIND=shum_real64), INTENT(OUT) :: weight_t_r(points) ! Weight applied to value at top ! right corner of source gridbox -REAL(KIND=real64), INTENT(OUT) :: weight_b_l(points) +REAL(KIND=shum_real64), INTENT(OUT) :: weight_b_l(points) ! Weight applied to value at bot. ! left corner of source gridbox -REAL(KIND=real64), INTENT(OUT) :: weight_b_r(points) +REAL(KIND=shum_real64), INTENT(OUT) :: weight_b_r(points) ! Weight applied to value at bot. ! right corner of source gridbox -REAL(KIND=real64), INTENT(OUT) :: weight_t_l(points) +REAL(KIND=shum_real64), INTENT(OUT) :: weight_t_l(points) ! Weight applied to value at top ! left corner of source gridbox -LOGICAL(KIND=bool), INTENT(IN) :: cyclic ! =T, then source data is cyclic - ! =F, then source data is non-cyclic +LOGICAL(KIND=shum_bool), INTENT(IN) :: cyclic ! =T, then source data is cyclic + ! =F, then source data is + ! non-cyclic !--- Local variables:--------------------------------------------------- -REAL(KIND=real64) :: t_lambda(points) ! Local value of target +REAL(KIND=shum_real64) :: t_lambda(points) ! Local value of target ! longitude -INTEGER(KIND=int64) :: ixp1(points) ! Longitudinal index plus 1 -INTEGER(KIND=int64) :: ix(points) ! Longitudinal index -INTEGER(KIND=int64) :: iy(points) ! Latitudinal index +INTEGER(KIND=shum_int64) :: ixp1(points) ! Longitudinal index plus 1 +INTEGER(KIND=shum_int64) :: ix(points) ! Longitudinal index +INTEGER(KIND=shum_int64) :: iy(points) ! Latitudinal index ! ---------------------------------------------------------------------- CALL f_shum_find_source_box_indices & @@ -236,7 +313,7 @@ SUBROUTINE f_shum_horizontal_field_bi_lin_interp_get_coeffs & , points_lambda_srce, points_phi_srce, points, cyclic & , t_lambda, ixp1, ix, iy ) -CALL f_shum_calc_weights & +CALL f_shum_calc_weights & ( weight_t_r, weight_b_r, weight_t_l, weight_b_l & , lambda_srce, phi_srce, phi_targ & , points_lambda_srce, points_phi_srce, points & @@ -254,54 +331,55 @@ SUBROUTINE f_shum_find_source_box_indices & IMPLICIT NONE -INTEGER(KIND=int64), INTENT(IN) :: points_lambda_srce +INTEGER(KIND=shum_int64), INTENT(IN) :: points_lambda_srce ! Number of lambda points on source grid -INTEGER(KIND=int64), INTENT(IN) :: points_phi_srce +INTEGER(KIND=shum_int64), INTENT(IN) :: points_phi_srce ! Number of phi points on source grid -INTEGER(KIND=int64), INTENT(IN) :: points +INTEGER(KIND=shum_int64), INTENT(IN) :: points ! Total number of points on target grid -INTEGER(KIND=int64), INTENT(OUT) :: index_b_l(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: index_b_l(points) ! Index of bottom left corner ! of source gridbox -INTEGER(KIND=int64), INTENT(OUT) :: index_b_r(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: index_b_r(points) ! Index of bottom right corner ! of source gridbox -INTEGER(KIND=int64), INTENT(OUT) :: index_t_l(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: index_t_l(points) ! Index of top left corner ! of source gridbox -INTEGER(KIND=int64), INTENT(OUT) :: index_t_r(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: index_t_r(points) ! Index of top right corner ! of source gridbox -INTEGER(KIND=int64), INTENT(OUT) :: ixp1(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: ixp1(points) ! Longitudinal index plus 1 -INTEGER(KIND=int64), INTENT(OUT) :: ix(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: ix(points) ! Longitudinal index -INTEGER(KIND=int64), INTENT(OUT) :: iy(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: iy(points) ! Latitudinal index -REAL(KIND=real64), INTENT(IN) :: lambda_targ(points) +REAL(KIND=shum_real64), INTENT(IN) :: lambda_targ(points) ! Lambda coords of target grid in degrees ! using same rotation as source grid -REAL(KIND=real64), INTENT(IN) :: phi_targ(points) +REAL(KIND=shum_real64), INTENT(IN) :: phi_targ(points) ! Phi coords of target grid in degrees using ! same rotation as source grid -REAL(KIND=real64), INTENT(IN) :: lambda_srce(points_lambda_srce) +REAL(KIND=shum_real64), INTENT(IN) :: lambda_srce(points_lambda_srce) ! Lambda coords of source grid in degrees -REAL(KIND=real64), INTENT(IN) :: phi_srce(points_phi_srce) +REAL(KIND=shum_real64), INTENT(IN) :: phi_srce(points_phi_srce) ! Phi coords of source grid in degrees -REAL(KIND=real64), INTENT(OUT) :: t_lambda(points) +REAL(KIND=shum_real64), INTENT(OUT) :: t_lambda(points) ! Local value of target longitude -LOGICAL(KIND=bool), INTENT(IN) :: cyclic ! =T, then source data is cyclic - ! =F, then source data is non-cyclic +LOGICAL(KIND=shum_bool), INTENT(IN) :: cyclic ! =T, then source data is cyclic + ! =F, then source data is + ! non-cyclic !--- Local variables:--------------------------------------------------- -INTEGER(KIND=int64) :: i ! Loop index +INTEGER(KIND=shum_int64) :: i ! Loop index ! Variables for divide-and-conquer search of latitude/longitude arrays -INTEGER(KIND=int64) :: iupr ! Uppermost-point -INTEGER(KIND=int64) :: ilwr ! Lowermost-point -INTEGER(KIND=int64) :: imid ! Mid-point +INTEGER(KIND=shum_int64) :: iupr ! Uppermost-point +INTEGER(KIND=shum_int64) :: ilwr ! Lowermost-point +INTEGER(KIND=shum_int64) :: imid ! Mid-point ! ---------------------------------------------------------------------- ! 1. Initialise arrays @@ -309,19 +387,19 @@ SUBROUTINE f_shum_find_source_box_indices & ! 1.1 Scale target longitude so that it falls between ! lambda_srce(1) and lambda_srce(1) + 360 DO i=1, points - t_lambda(i) = MOD(((lambda_targ(i)-lambda_srce(1))+720.0_real64), & - 360.0_real64) + lambda_srce(1) + t_lambda(i) = MOD(((lambda_targ(i)-lambda_srce(1))+720.0_shum_real64), & + 360.0_shum_real64) + lambda_srce(1) END DO IF (cyclic) THEN - DO i=1_int64, points - ix(i) = 0_int64 - iy(i) = 1_int64 + DO i=one64, points + ix(i) = zero64 + iy(i) = one64 END DO ELSE - DO i=1_int64, points - ix(i) = 1_int64 - iy(i) = 1_int64 + DO i=one64, points + ix(i) = one64 + iy(i) = one64 END DO END IF @@ -331,18 +409,18 @@ SUBROUTINE f_shum_find_source_box_indices & ! Longitude !$OMP PARALLEL DO SCHEDULE(STATIC) & -!$OMP& SHARED(points,points_lambda_srce,t_lambda,lambda_srce) & -!$OMP& SHARED(points_phi_srce,phi_srce,phi_targ,ix,iy) & -!$OMP& PRIVATE(i,iupr,ilwr,imid) DEFAULT(NONE) -DO i=1_int64, points +!$OMP SHARED(points,points_lambda_srce,t_lambda,lambda_srce) & +!$OMP SHARED(points_phi_srce,phi_srce,phi_targ,ix,iy) & +!$OMP PRIVATE(i,iupr,ilwr,imid) DEFAULT(NONE) +DO i=one64, points ! Divide and conquer should more efficient than a brute-force loop over ! all i - ilwr = 0_int64 !first point(-1 in case we are cyclic) - iupr = points_lambda_srce + 1_int64 !last point (+1 in case we are cyclic) - DO WHILE ( (iupr-ilwr) > 1_int64 ) - imid = (ilwr+iupr)/2_int64 + ilwr = zero64 ! first point (-1 in case we are cyclic) + iupr = points_lambda_srce + one64 ! last point (+1 in case we are cyclic) + DO WHILE ( (iupr-ilwr) > one64 ) + imid = (ilwr+iupr)/two64 IF ( lambda_srce(imid) > t_lambda(i) ) THEN iupr = imid END IF @@ -353,10 +431,10 @@ SUBROUTINE f_shum_find_source_box_indices & ix(i) = ilwr - ilwr = 1_int64 !first point + ilwr = one64 !first point iupr = points_phi_srce !last point - DO WHILE ( (iupr-ilwr) > 1_int64 ) - imid = (ilwr+iupr)/2_int64 + DO WHILE ( (iupr-ilwr) > one64 ) + imid = (ilwr+iupr)/two64 IF ( phi_srce(imid) > phi_targ(i) ) THEN iupr = imid END IF @@ -385,20 +463,20 @@ SUBROUTINE f_shum_find_source_box_indices & ! but we shall keep it here in case). IF (ix(i) < 1) THEN ix(i) = points_lambda_srce - t_lambda(i) = t_lambda(i) + 360.0_real64 + t_lambda(i) = t_lambda(i) + 360.0_shum_real64 END IF ! Set index for one sided difference if target point to north or ! south of source area. - iy(i) = MAX(iy(i),1_int64) - iy(i) = MIN(iy(i),points_phi_srce-1_int64) + iy(i) = MAX(iy(i),one64) + iy(i) = MIN(iy(i),points_phi_srce-one64) ! 2-D indices - index_b_l(i) = ix(i)+(iy(i)-1_int64)*points_lambda_srce - index_b_r(i) = index_b_l(i)+1_int64 + index_b_l(i) = ix(i)+(iy(i)-one64)*points_lambda_srce + index_b_r(i) = index_b_l(i)+one64 ! Correct for cyclic boundaries if target point outside source grid. - ixp1(i) = ix(i) + 1_int64 + ixp1(i) = ix(i) + one64 IF (ix(i) == points_lambda_srce) THEN index_b_r(i) = index_b_r(i) - points_lambda_srce ixp1(i) = ixp1(i) - points_lambda_srce @@ -424,33 +502,33 @@ SUBROUTINE f_shum_find_source_box_indices & ! less than lambda_srce(1). I can't think of a legitimate case where this ! would be the case and it needs further investigating. IF (ix(i) == points_lambda_srce) THEN - IF (ABS(t_lambda(i)-lambda_srce(1)-360.0_real64) < & + IF (ABS(t_lambda(i)-lambda_srce(1)-360.0_shum_real64) < & t_lambda(i)-lambda_srce(ix(i))) THEN - t_lambda(i) = t_lambda(i) - 360.0_real64 + t_lambda(i) = t_lambda(i) - 360.0_shum_real64 ix(i) = 1 END IF END IF ! Set index for one sided difference if outside source area - ix(i) = MAX(ix(i),1_int64) - ix(i) = MIN(ix(i),points_lambda_srce-1_int64) - IF (ix(i) < 1_int64) THEN ! IX(I) < 1 if POINTS_LAMBDA_SRCE = 1 - ix(i) = 1_int64 + ix(i) = MAX(ix(i),one64) + ix(i) = MIN(ix(i),points_lambda_srce-one64) + IF (ix(i) < one64) THEN ! IX(I) < 1 if POINTS_LAMBDA_SRCE = 1 + ix(i) = one64 END IF - ixp1(i) = ix(i) + 1_int64 + ixp1(i) = ix(i) + one64 ixp1(i) = MIN(ixp1(i),points_lambda_srce) ! Set index for one sided difference if outside source area - iy(i) = MAX(iy(i),1_int64) - iy(i) = MIN(iy(i),points_phi_srce-1_int64) - IF (iy(i) < 1_int64) THEN ! IY(I) < 1 if POINTS_PHI_SRCE = 1 - iy(i) = 1_int64 + iy(i) = MAX(iy(i),one64) + iy(i) = MIN(iy(i),points_phi_srce-one64) + IF (iy(i) < one64) THEN ! IY(I) < 1 if POINTS_PHI_SRCE = 1 + iy(i) = one64 END IF ! 2-D indices - index_b_l(i) = ix(i) + (iy(i)-1_int64)*points_lambda_srce - index_b_r(i) = ixp1(i) + (iy(i)-1_int64)*points_lambda_srce + index_b_l(i) = ix(i) + (iy(i)-one64)*points_lambda_srce + index_b_r(i) = ixp1(i) + (iy(i)-one64)*points_lambda_srce index_t_l(i) = index_b_l(i) + points_lambda_srce index_t_r(i) = index_b_r(i) + points_lambda_srce @@ -470,72 +548,73 @@ SUBROUTINE f_shum_calc_weights & IMPLICIT NONE -INTEGER(KIND=int64), INTENT(IN) :: points_lambda_srce +INTEGER(KIND=shum_int64), INTENT(IN) :: points_lambda_srce ! Number of lambda points on source grid -INTEGER(KIND=int64), INTENT(IN) :: points_phi_srce +INTEGER(KIND=shum_int64), INTENT(IN) :: points_phi_srce ! Number of phi points on source grid -INTEGER(KIND=int64), INTENT(IN) :: points +INTEGER(KIND=shum_int64), INTENT(IN) :: points ! Total number of points on target grid -INTEGER(KIND=int64), INTENT(IN) :: ixp1(points) +INTEGER(KIND=shum_int64), INTENT(IN) :: ixp1(points) ! Longitudinal index plus 1 -INTEGER(KIND=int64), INTENT(IN) :: ix(points) +INTEGER(KIND=shum_int64), INTENT(IN) :: ix(points) ! Longitudinal index -INTEGER(KIND=int64), INTENT(IN) :: iy(points) +INTEGER(KIND=shum_int64), INTENT(IN) :: iy(points) ! Latitudinal index -REAL(KIND=real64), INTENT(IN) :: phi_targ(points) +REAL(KIND=shum_real64), INTENT(IN) :: phi_targ(points) ! Phi coords of target grid in degrees using ! same rotation as source grid -REAL(KIND=real64), INTENT(IN) :: lambda_srce(points_lambda_srce) +REAL(KIND=shum_real64), INTENT(IN) :: lambda_srce(points_lambda_srce) ! Lambda coords of source grid in degrees -REAL(KIND=real64), INTENT(IN) :: phi_srce(points_phi_srce) +REAL(KIND=shum_real64), INTENT(IN) :: phi_srce(points_phi_srce) ! Phi coords of source grid in degrees -REAL(KIND=real64), INTENT(IN) :: t_lambda(points) +REAL(KIND=shum_real64), INTENT(IN) :: t_lambda(points) ! Local value of target longitude -REAL(KIND=real64), INTENT(OUT) :: weight_t_r(points) +REAL(KIND=shum_real64), INTENT(OUT) :: weight_t_r(points) ! Weight applied to value at top ! right corner of source gridbox -REAL(KIND=real64), INTENT(OUT) :: weight_b_l(points) +REAL(KIND=shum_real64), INTENT(OUT) :: weight_b_l(points) ! Weight applied to value at bot. ! left corner of source gridbox -REAL(KIND=real64), INTENT(OUT) :: weight_b_r(points) +REAL(KIND=shum_real64), INTENT(OUT) :: weight_b_r(points) ! Weight applied to value at bot. ! right corner of source gridbox -REAL(KIND=real64), INTENT(OUT) :: weight_t_l(points) +REAL(KIND=shum_real64), INTENT(OUT) :: weight_t_l(points) ! Weight applied to value at top ! left corner of source gridbox !--- Local variables:--------------------------------------------------- -REAL(KIND=real64) :: a ! Longitudinal weight -REAL(KIND=real64) :: b ! Latitudinal weight -INTEGER(KIND=int64) :: i ! Loop index +REAL(KIND=shum_real64) :: a ! Longitudinal weight +REAL(KIND=shum_real64) :: b ! Latitudinal weight +INTEGER(KIND=shum_int64) :: i ! Loop index ! ---------------------------------------------------------------------- ! 1. Compute interpolation weights DO i=1, points ! Calculate basic weights - a = (MOD(360.0_real64+lambda_srce(ixp1(i))-lambda_srce(ix(i)),360.0_real64)) - IF (a /= 0.0_real64) THEN + a = MOD(360.0_shum_real64 + lambda_srce(ixp1(i)) - lambda_srce(ix(i)), & + 360.0_shum_real64) + IF (a /= 0.0_shum_real64) THEN ! If t_lambda - lambda_source is negative then just copy last value across. - a = (MAX(t_lambda(i)-lambda_srce(ix(i)),0.0_real64))/a + a = (MAX(t_lambda(i)-lambda_srce(ix(i)),0.0_shum_real64))/a ELSE - a = 0.0_real64 + a = 0.0_shum_real64 END IF ! If we only have 1 row then we need to make sure we can cope. b = ABS(phi_srce(iy(i))-phi_srce(MIN(iy(i)+1,points_phi_srce))) - IF (b /= 0.0_real64) THEN - b = MAX(phi_targ(i)-phi_srce(iy(i)),0.0_real64)/b + IF (b /= 0.0_shum_real64) THEN + b = MAX(phi_targ(i)-phi_srce(iy(i)),0.0_shum_real64)/b ELSE - b = 0.0_real64 + b = 0.0_shum_real64 END IF ! Calculate bi-linear interpolation weights weight_t_r(i) = a*b - weight_b_l(i) = (1.0_real64-a)*(1.0_real64-b) - weight_t_l(i) = (1.0_real64-a)*b - weight_b_r(i) = a*(1.0_real64-b) + weight_b_l(i) = (1.0_shum_real64-a)*(1.0_shum_real64-b) + weight_t_l(i) = (1.0_shum_real64-a)*b + weight_b_r(i) = a*(1.0_shum_real64-b) END DO @@ -577,62 +656,62 @@ SUBROUTINE f_shum_cart_horizontal_field_bi_lin_interp_get_coeffs & IMPLICIT NONE -INTEGER(KIND=int64), INTENT(IN) :: points_x_srce +INTEGER(KIND=shum_int64), INTENT(IN) :: points_x_srce ! Number of X on source grid -INTEGER(KIND=int64), INTENT(IN) :: points_y_srce +INTEGER(KIND=shum_int64), INTENT(IN) :: points_y_srce ! Number of Y points on source grid -INTEGER(KIND=int64), INTENT(IN) :: points +INTEGER(KIND=shum_int64), INTENT(IN) :: points ! Total number of points on target grid -INTEGER(KIND=int64), INTENT(IN) :: grid_type +INTEGER(KIND=shum_int64), INTENT(IN) :: grid_type ! Source grid type ! 3 - LAM no wray around ! 4 - LAM bicylic wrap in X and Y ! - LAM channel wrap in X only -INTEGER(KIND=int64), INTENT(OUT) :: index_b_l(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: index_b_l(points) ! Index of bottom left corner ! of source gridbox -INTEGER(KIND=int64), INTENT(OUT) :: index_b_r(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: index_b_r(points) ! Index of bottom right corner ! of source gridbox -INTEGER(KIND=int64), INTENT(OUT) :: index_t_l(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: index_t_l(points) ! Index of top left corner ! of source gridbox -INTEGER(KIND=int64), INTENT(OUT) :: index_t_r(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: index_t_r(points) ! Index of top right corner ! of source gridbox -REAL(KIND=real64), INTENT(IN) :: x_targ(points) +REAL(KIND=shum_real64), INTENT(IN) :: x_targ(points) ! X coords of target grid in m -REAL(KIND=real64), INTENT(IN) :: y_targ(points) +REAL(KIND=shum_real64), INTENT(IN) :: y_targ(points) ! Y coords of target grid in m -REAL(KIND=real64), INTENT(IN) :: x_srce(points_x_srce) +REAL(KIND=shum_real64), INTENT(IN) :: x_srce(points_x_srce) ! X coords of source grid in m -REAL(KIND=real64), INTENT(IN) :: y_srce(points_y_srce) +REAL(KIND=shum_real64), INTENT(IN) :: y_srce(points_y_srce) ! Y coords of source grid in m -REAL(KIND=real64), INTENT(OUT) :: weight_t_r(points) +REAL(KIND=shum_real64), INTENT(OUT) :: weight_t_r(points) ! Weight applied to value at top ! right corner of source gridbox -REAL(KIND=real64), INTENT(OUT) :: weight_b_l(points) +REAL(KIND=shum_real64), INTENT(OUT) :: weight_b_l(points) ! Weight applied to value at bot. ! left corner of source gridbox -REAL(KIND=real64), INTENT(OUT) :: weight_b_r(points) +REAL(KIND=shum_real64), INTENT(OUT) :: weight_b_r(points) ! Weight applied to value at bot. ! right corner of source gridbox -REAL(KIND=real64), INTENT(OUT) :: weight_t_l(points) +REAL(KIND=shum_real64), INTENT(OUT) :: weight_t_l(points) ! Weight applied to value at top ! left corner of source gridbox -INTEGER(KIND=int64), INTENT(OUT) :: icode ! error code +INTEGER(KIND=shum_int64), INTENT(OUT) :: icode ! error code CHARACTER(LEN=*), INTENT(OUT) :: cmessage ! error message !--- Local variables:--------------------------------------------------- -INTEGER(KIND=int64) :: ixp1(points) ! X index plus 1 -INTEGER(KIND=int64) :: ix(points) ! X index -INTEGER(KIND=int64) :: iyp1(points) ! Y index plus 1 -INTEGER(KIND=int64) :: iy(points) ! Y index +INTEGER(KIND=shum_int64) :: ixp1(points) ! X index plus 1 +INTEGER(KIND=shum_int64) :: ix(points) ! X index +INTEGER(KIND=shum_int64) :: iyp1(points) ! Y index plus 1 +INTEGER(KIND=shum_int64) :: iy(points) ! Y index ! ---------------------------------------------------------------------- CALL f_shum_find_source_cart_box_indices & @@ -663,72 +742,72 @@ SUBROUTINE f_shum_find_source_cart_box_indices & IMPLICIT NONE -INTEGER(KIND=int64), INTENT(IN) :: points_x_srce +INTEGER(KIND=shum_int64), INTENT(IN) :: points_x_srce ! Number of x points on source grid -INTEGER(KIND=int64), INTENT(IN) :: points_y_srce +INTEGER(KIND=shum_int64), INTENT(IN) :: points_y_srce ! Number of y points on source grid -INTEGER(KIND=int64), INTENT(IN) :: points +INTEGER(KIND=shum_int64), INTENT(IN) :: points ! Total number of points on target grid -INTEGER(KIND=int64), INTENT(IN) :: grid_type +INTEGER(KIND=shum_int64), INTENT(IN) :: grid_type ! Source grid type ! 3 - LAM no wray around ! 4 - LAM bicylic wrap in X and Y ! - LAM channel wrap in X only -INTEGER(KIND=int64), INTENT(OUT) :: index_b_l(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: index_b_l(points) ! Index of bottom left corner ! of source gridbox -INTEGER(KIND=int64), INTENT(OUT) :: index_b_r(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: index_b_r(points) ! Index of bottom right corner ! of source gridbox -INTEGER(KIND=int64), INTENT(OUT) :: index_t_l(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: index_t_l(points) ! Index of top left corner ! of source gridbox -INTEGER(KIND=int64), INTENT(OUT) :: index_t_r(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: index_t_r(points) ! Index of top right corner ! of source gridbox -INTEGER(KIND=int64), INTENT(OUT) :: ixp1(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: ixp1(points) ! X index plus 1 -INTEGER(KIND=int64), INTENT(OUT) :: ix(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: ix(points) ! X index -INTEGER(KIND=int64), INTENT(OUT) :: iyp1(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: iyp1(points) ! Y index plus 1 -INTEGER(KIND=int64), INTENT(OUT) :: iy(points) +INTEGER(KIND=shum_int64), INTENT(OUT) :: iy(points) ! Y index -INTEGER(KIND=int64), INTENT(OUT) :: icode ! error code +INTEGER(KIND=shum_int64), INTENT(OUT) :: icode ! error code CHARACTER(LEN=*), INTENT(OUT) :: cmessage ! error message -REAL(KIND=real64), INTENT(IN) :: x_targ(points) +REAL(KIND=shum_real64), INTENT(IN) :: x_targ(points) ! x coords of target grid in m -REAL(KIND=real64), INTENT(IN) :: y_targ(points) +REAL(KIND=shum_real64), INTENT(IN) :: y_targ(points) ! y coords of target grid in m -REAL(KIND=real64), INTENT(IN) :: x_srce(points_x_srce) +REAL(KIND=shum_real64), INTENT(IN) :: x_srce(points_x_srce) ! x coords of source grid in m -REAL(KIND=real64), INTENT(IN) :: y_srce(points_y_srce) +REAL(KIND=shum_real64), INTENT(IN) :: y_srce(points_y_srce) ! y coords of source grid in m !--- Local variables:--------------------------------------------------- -INTEGER(KIND=int64) :: i ! Loop index -INTEGER(KIND=int64) :: j ! Loop index +INTEGER(KIND=shum_int64) :: i ! Loop index +INTEGER(KIND=shum_int64) :: j ! Loop index ! Grid types -INTEGER(KIND=int64), PARAMETER :: bicylic = 4 -INTEGER(KIND=int64), PARAMETER :: channel = 3 -INTEGER(KIND=int64), PARAMETER :: nowrap = 2 +INTEGER(KIND=shum_int64), PARAMETER :: bicylic = 4 +INTEGER(KIND=shum_int64), PARAMETER :: channel = 3 +INTEGER(KIND=shum_int64), PARAMETER :: nowrap = 2 ! ---------------------------------------------------------------------- ! Set error code to zero -icode = 0_int64 +icode = zero64 ! 1. Initialise arrays IF (grid_type == bicylic) THEN - DO i = 1_int64, points - DO j = 1_int64, points_x_srce-1_int64 - IF ( x_targ(i) >= x_srce(j) .AND. x_targ(i) < x_srce(j+1_int64) ) THEN + DO i = one64, points + DO j = one64, points_x_srce-one64 + IF ( x_targ(i) >= x_srce(j) .AND. x_targ(i) < x_srce(j+one64) ) THEN ix(i) = j - ixp1(i) = j+1_int64 + ixp1(i) = j+one64 END IF END DO ! Wrap around grid @@ -743,51 +822,51 @@ SUBROUTINE f_shum_find_source_cart_box_indices & ! 0.0 max_x_srce ! ! Test for points in locations A or B - IF ( x_targ(i) < x_srce(1_int64) .OR. & + IF ( x_targ(i) < x_srce(one64) .OR. & x_targ(i) >= x_srce(points_x_srce) ) THEN ix(i) = points_x_srce - ixp1(i) = 1_int64 + ixp1(i) = one64 END IF - DO j = 1_int64, points_y_srce-1_int64 + DO j = one64, points_y_srce-one64 IF ( y_targ(i) >= y_srce(j) .AND. y_targ(i) < y_srce(j+1) ) THEN iy(i) = j - iyp1(i) = j+1_int64 + iyp1(i) = j+one64 END IF END DO ! Wrap around tests as for X - IF ( y_targ(i) < y_srce(1_int64) .OR. & + IF ( y_targ(i) < y_srce(one64) .OR. & y_targ(i) >= y_srce(points_y_srce) ) THEN iy(i) = points_y_srce - iyp1(i) = 1_int64 + iyp1(i) = one64 END IF END DO ELSE IF (grid_type == nowrap) THEN - DO i = 1_int64, points - DO j = 1_int64, points_x_srce-1_int64 - IF ( x_targ(i) >= x_srce(j) .AND. x_targ(i) < x_srce(j+1_int64) ) THEN + DO i = one64, points + DO j = one64, points_x_srce-one64 + IF ( x_targ(i) >= x_srce(j) .AND. x_targ(i) < x_srce(j+one64) ) THEN ix(i) = j - ixp1(i) = j+1_int64 + ixp1(i) = j+one64 END IF END DO - IF ( x_targ(i) < x_srce(1_int64) ) THEN - ix(i) = 1_int64 - ixp1(i) = 1_int64 + IF ( x_targ(i) < x_srce(one64) ) THEN + ix(i) = one64 + ixp1(i) = one64 END IF IF ( x_targ(i) >= x_srce(points_x_srce) ) THEN ix(i) = points_x_srce ixp1(i) = points_x_srce END IF - DO j = 1_int64, points_y_srce-1 - IF ( y_targ(i) >= y_srce(j) .AND. y_targ(i) < y_srce(j+1_int64) ) THEN + DO j = one64, points_y_srce-1 + IF ( y_targ(i) >= y_srce(j) .AND. y_targ(i) < y_srce(j+one64) ) THEN iy(i) = j - iyp1(i) = j+1_int64 + iyp1(i) = j+one64 END IF END DO - IF ( y_targ(i) < y_srce(1_int64) ) THEN - iy(i) = 1_int64 - iyp1(i) = 1_int64 + IF ( y_targ(i) < y_srce(one64) ) THEN + iy(i) = one64 + iyp1(i) = one64 END IF IF ( y_targ(i) >= y_srce(points_y_srce) ) THEN iy(i) = points_y_srce @@ -797,29 +876,29 @@ SUBROUTINE f_shum_find_source_cart_box_indices & ELSE IF (grid_type == channel) THEN ! wrap in x direction only - DO i = 1_int64, points - DO j = 1_int64, points_x_srce-1_int64 - IF ( x_targ(i) >= x_srce(j) .AND. x_targ(i) < x_srce(j+1_int64) ) THEN + DO i = one64, points + DO j = one64, points_x_srce-one64 + IF ( x_targ(i) >= x_srce(j) .AND. x_targ(i) < x_srce(j+one64) ) THEN ix(i) = j - ixp1(i) = j+1_int64 + ixp1(i) = j+one64 END IF END DO ! Wrap around tests as for bicyclic case IF ( x_targ(i) < x_srce(1) .OR. x_targ(i) >= x_srce(points_x_srce) ) THEN ix(i) = points_x_srce - ixp1(i) = 1_int64 + ixp1(i) = one64 END IF ! No wrap around - DO j = 1_int64, points_y_srce-1 - IF ( y_targ(i) >= y_srce(j) .AND. y_targ(i) < y_srce(j+1_int64) ) THEN + DO j = one64, points_y_srce-1 + IF ( y_targ(i) >= y_srce(j) .AND. y_targ(i) < y_srce(j+one64) ) THEN iy(i) = j - iyp1(i) = j+1_int64 + iyp1(i) = j+one64 END IF END DO - IF (y_targ(i) < y_srce(1_int64) ) THEN - iy(i) = 1_int64 - iyp1(i) = 1_int64 + IF (y_targ(i) < y_srce(one64) ) THEN + iy(i) = one64 + iyp1(i) = one64 END IF IF (y_targ(i) >= y_srce(points_y_srce) ) THEN iy(i) = points_y_srce @@ -847,11 +926,11 @@ SUBROUTINE f_shum_find_source_cart_box_indices & ! Right hand bottom corner coordinates ixp1, iy ! Right hand top corner coordinates ixp1, iyp1 -DO i = 1_int64, points - index_b_l(i) = ix(i) + (iy(i)-1_int64)*points_x_srce - index_b_r(i) = ixp1(i) + (iy(i)-1_int64)*points_x_srce - index_t_l(i) = ix(i) + (iyp1(i)-1_int64)*points_x_srce - index_t_r(i) = ixp1(i) + (iyp1(i)-1_int64)*points_x_srce +DO i = one64, points + index_b_l(i) = ix(i) + (iy(i)-one64)*points_x_srce + index_b_r(i) = ixp1(i) + (iy(i)-one64)*points_x_srce + index_t_l(i) = ix(i) + (iyp1(i)-one64)*points_x_srce + index_t_r(i) = ixp1(i) + (iyp1(i)-one64)*points_x_srce END DO @@ -868,67 +947,67 @@ SUBROUTINE f_shum_calc_cart_weights & IMPLICIT NONE -INTEGER(KIND=int64), INTENT(IN) :: points_x_srce +INTEGER(KIND=shum_int64), INTENT(IN) :: points_x_srce ! Number of x points on source grid -INTEGER(KIND=int64), INTENT(IN) :: points_y_srce +INTEGER(KIND=shum_int64), INTENT(IN) :: points_y_srce ! Number of y points on source grid -INTEGER(KIND=int64), INTENT(IN) :: points +INTEGER(KIND=shum_int64), INTENT(IN) :: points ! Total number of points on target grid -INTEGER(KIND=int64), INTENT(IN) :: ixp1(points) +INTEGER(KIND=shum_int64), INTENT(IN) :: ixp1(points) ! X index plus 1 -INTEGER(KIND=int64), INTENT(IN) :: ix(points) +INTEGER(KIND=shum_int64), INTENT(IN) :: ix(points) ! X index -INTEGER(KIND=int64), INTENT(IN) :: iyp1(points) +INTEGER(KIND=shum_int64), INTENT(IN) :: iyp1(points) ! Y index plus 1 -INTEGER(KIND=int64), INTENT(IN) :: iy(points) +INTEGER(KIND=shum_int64), INTENT(IN) :: iy(points) ! Y index -REAL(KIND=real64), INTENT(IN) :: x_targ(points) +REAL(KIND=shum_real64), INTENT(IN) :: x_targ(points) ! X coords of target grid in m -REAL(KIND=real64), INTENT(IN) :: y_targ(points) +REAL(KIND=shum_real64), INTENT(IN) :: y_targ(points) ! y coords of target grid in m -REAL(KIND=real64), INTENT(IN) :: x_srce(points_x_srce) +REAL(KIND=shum_real64), INTENT(IN) :: x_srce(points_x_srce) ! x coords of source grid in m -REAL(KIND=real64), INTENT(IN) :: y_srce(points_y_srce) +REAL(KIND=shum_real64), INTENT(IN) :: y_srce(points_y_srce) ! y coords of source grid in m -REAL(KIND=real64), INTENT(OUT) :: weight_t_r(points) +REAL(KIND=shum_real64), INTENT(OUT) :: weight_t_r(points) ! Weight applied to value at top ! right corner of source gridbox -REAL(KIND=real64), INTENT(OUT) :: weight_b_l(points) +REAL(KIND=shum_real64), INTENT(OUT) :: weight_b_l(points) ! Weight applied to value at bot. ! left corner of source gridbox -REAL(KIND=real64), INTENT(OUT) :: weight_b_r(points) +REAL(KIND=shum_real64), INTENT(OUT) :: weight_b_r(points) ! Weight applied to value at bot. ! right corner of source gridbox -REAL(KIND=real64), INTENT(OUT) :: weight_t_l(points) +REAL(KIND=shum_real64), INTENT(OUT) :: weight_t_l(points) ! Weight applied to value at top ! left corner of source gridbox !--- Local variables:--------------------------------------------------- -REAL(KIND=real64) :: a ! x weight -REAL(KIND=real64) :: b ! y weight -REAL(KIND=real64) :: dx_srce ! x grid length on source grid (m) -REAL(KIND=real64) :: dy_srce ! y grid length on source grid (m) -REAL(KIND=real64) :: max_x_srce ! x length of source grid (m) -REAL(KIND=real64) :: max_y_srce ! y length of source grid (m) -INTEGER(KIND=int64) :: i ! Loop index +REAL(KIND=shum_real64) :: a ! x weight +REAL(KIND=shum_real64) :: b ! y weight +REAL(KIND=shum_real64) :: dx_srce ! x grid length on source grid (m) +REAL(KIND=shum_real64) :: dy_srce ! y grid length on source grid (m) +REAL(KIND=shum_real64) :: max_x_srce ! x length of source grid (m) +REAL(KIND=shum_real64) :: max_y_srce ! y length of source grid (m) +INTEGER(KIND=shum_int64) :: i ! Loop index ! ---------------------------------------------------------------------- ! Cartesian grid spacing assumed regular -dx_srce = x_srce(2_int64) - x_srce(1_int64) -dy_srce = y_srce(2_int64) - y_srce(1_int64) +dx_srce = x_srce(two64) - x_srce(one64) +dy_srce = y_srce(two64) - y_srce(one64) ! Required for wrap around cases max_x_srce = points_x_srce*dx_srce max_y_srce = points_y_srce*dy_srce ! 1. Compute interpolation weights -DO i=1_int64, points +DO i=one64, points ! Calculate basic weights need different ix values IF (ixp1(i) > ix(i)) THEN - a = (MAX(x_targ(i)-x_srce(ix(i)),0.0_real64))/dx_srce + a = (MAX(x_targ(i)-x_srce(ix(i)),0.0_shum_real64))/dx_srce ELSE IF (ixp1(i) < ix(i)) THEN ! Wrap around case ! Note first grid point X is not always 0.0 depends on whether U, V or P ! grid so target grid value can be > final source_grid_y but less than @@ -941,33 +1020,33 @@ SUBROUTINE f_shum_calc_cart_weights & ! The new point in the wrap around region can be in position A or B IF ( x_targ(i) < x_srce(ix(i))) THEN ! target in position A - a = (MAX(x_targ(i)-(x_srce(ix(i))-max_x_srce),0.0_real64))/dx_srce + a = (MAX(x_targ(i)-(x_srce(ix(i))-max_x_srce),0.0_shum_real64))/dx_srce ELSE ! Target in position B - a = (MAX(x_targ(i)-x_srce(ix(i)),0.0_real64))/dx_srce + a = (MAX(x_targ(i)-x_srce(ix(i)),0.0_shum_real64))/dx_srce END IF ELSE ! values equal so zero - a = 0.0_real64 + a = 0.0_shum_real64 END IF ! Calculate basic weights need different iy values IF (iyp1(i) > iy(i)) THEN - b = (MAX(y_targ(i)-y_srce(iy(i)),0.0_real64))/dy_srce + b = (MAX(y_targ(i)-y_srce(iy(i)),0.0_shum_real64))/dy_srce ELSE IF (iyp1(i) < iy(i)) THEN ! Wrap around case ! Same problems in Y direction as X direction. IF ( y_targ(i) < y_srce(iy(i))) THEN ! Like position A - b = (MAX(y_targ(i)-(y_srce(iy(i))-max_y_srce),0.0_real64))/dy_srce + b = (MAX(y_targ(i)-(y_srce(iy(i))-max_y_srce),0.0_shum_real64))/dy_srce ELSE ! like position B - b = (MAX(y_targ(i)-y_srce(iy(i)),0.0_real64))/dy_srce + b = (MAX(y_targ(i)-y_srce(iy(i)),0.0_shum_real64))/dy_srce END IF ELSE ! Corners the same so zero - b = 0.0_real64 + b = 0.0_shum_real64 END IF ! Calculate bi-linear interpolation weights weight_t_r(i) = a*b - weight_b_l(i) = (1.0_real64-a)*(1.0_real64-b) - weight_t_l(i) = (1.0_real64-a)*b - weight_b_r(i) = a*(1.0_real64-b) + weight_b_l(i) = (1.0_shum_real64-a)*(1.0_shum_real64-b) + weight_t_l(i) = (1.0_shum_real64-a)*b + weight_b_r(i) = a*(1.0_shum_real64-b) END DO diff --git a/shum_spiral_search/src/f_shum_spiral_search.f90 b/shum_spiral_search/src/f_shum_spiral_search.f90 index f9e870c..ecf094a 100644 --- a/shum_spiral_search/src/f_shum_spiral_search.f90 +++ b/shum_spiral_search/src/f_shum_spiral_search.f90 @@ -197,7 +197,6 @@ FUNCTION f_shum_spiral_arg64 & ! Loop over every unresolved point. DO k=1, no_point_unres - distance(:,:)=rMDI ! Find i/j index for point. unres_j = (index_unres(k)-1_int64)/points_lambda + 1_int64 unres_i = index_unres(k)-(unres_j-1_int64)*points_lambda @@ -239,8 +238,6 @@ FUNCTION f_shum_spiral_arg64 & curr_dist_valid_min = (planet_radius*shum_pi_const)+1.0_real64 ! We want nearest one so we dont want to limit search to max_dist. curr_dist_invalid_min = rMDI -! Give some initial rMDI values in case they are never set - distance(:,:)=rMDI found=.FALSE. @@ -260,58 +257,42 @@ FUNCTION f_shum_spiral_arg64 & east=-1_int64 ! Find how many points can go south and be inside search_dist - DO j = 1, points_phi - IF (unres_j-j > 0) THEN - tempdist=calc_distance(lats(unres_j),lons(unres_i), & - lats(unres_j-j),lons(unres_i), & - planet_radius) - IF (tempdist > search_dist) THEN - south=j-1_int64 - EXIT - END IF - ELSE + DO j = 1, unres_j-1 + tempdist=calc_distance(lats(unres_j),lons(unres_i), & + lats(unres_j-j),lons(unres_i), & + planet_radius) + IF (tempdist > search_dist) THEN + south=j-1_int64 EXIT END IF END DO ! Find how many points can go north and be inside search_dist - DO j = 1, points_phi - IF (unres_j+j <= points_phi) THEN - tempdist=calc_distance(lats(unres_j),lons(unres_i), & - lats(unres_j+j),lons(unres_i), & - planet_radius) - IF (tempdist > search_dist) THEN - north=j-1_int64 - EXIT - END IF - ELSE + DO j = 1, points_phi-unres_j + tempdist=calc_distance(lats(unres_j),lons(unres_i), & + lats(unres_j+j),lons(unres_i), & + planet_radius) + IF (tempdist > search_dist) THEN + north=j-1_int64 EXIT END IF END DO ! Find how many points can go west and be inside search_dist - DO i = 1, points_lambda - IF (unres_i-i > 0) THEN - tempdist=calc_distance(lats(unres_j),lons(unres_i), & - lats(unres_j),lons(unres_i-i), & - planet_radius) - IF (tempdist > search_dist) THEN - west=i-1_int64 - EXIT - END IF - ELSE + DO i = 1, unres_i-1 + tempdist=calc_distance(lats(unres_j),lons(unres_i), & + lats(unres_j),lons(unres_i-i), & + planet_radius) + IF (tempdist > search_dist) THEN + west=i-1_int64 EXIT END IF END DO ! Find how many points can go east and be inside search_dist - DO i = 1, points_lambda - IF (unres_i+i <= points_lambda) THEN - tempdist=calc_distance(lats(unres_j),lons(unres_i), & - lats(unres_j),lons(unres_i+i), & - planet_radius) - IF (tempdist > search_dist) THEN - east=i-1_int64 - EXIT - END IF - ELSE + DO i = 1, points_lambda-unres_i + tempdist=calc_distance(lats(unres_j),lons(unres_i), & + lats(unres_j),lons(unres_i+i), & + planet_radius) + IF (tempdist > search_dist) THEN + east=i-1_int64 EXIT END IF END DO @@ -343,6 +324,13 @@ FUNCTION f_shum_spiral_arg64 & ! Calculate distance from point distance(i,j) = calc_distance(lats(unres_j), lons(unres_i), & lats(j), lons(i), planet_radius) + END IF + END DO + + DO i = 1, points_lambda + l = i+(j - 1_int64)*points_lambda + ! Check to see if it is a resolved point + IF (.NOT. unres_mask(l)) THEN ! Same type of point. IF (lsm(l) .EQV. is_land_field) THEN ! If current distance is less than any previous min store it. @@ -405,6 +393,13 @@ FUNCTION f_shum_spiral_arg64 & ! Calculate distance from point distance(i,j) = calc_distance(lats(unres_j),lons(unres_i), & lats(j),lons(i), planet_radius) + END IF + END DO + + DO i = unres_i-west, unres_i+east + l = i+(j - 1_int64)*points_lambda + ! Check to see if it is a resolved point + IF (.NOT. unres_mask(l)) THEN ! Same type of point and distance less than maximum distance. IF ((lsm(l) .EQV. is_land_field) .AND. & (distance(i,j) < search_dist)) THEN @@ -603,7 +598,6 @@ FUNCTION f_shum_spiral_arg32 & ! Loop over every unresolved point. DO k=1, no_point_unres - distance(:,:)=rMDI_32b ! Find i/j index for point. unres_j = (index_unres(k)-1_int32)/points_lambda + 1_int32 unres_i = index_unres(k)-(unres_j-1_int32)*points_lambda @@ -645,8 +639,6 @@ FUNCTION f_shum_spiral_arg32 & curr_dist_valid_min = (planet_radius*shum_pi_const_32)+1.0_real32 ! We want nearest one so we dont want to limit search to max_dist. curr_dist_invalid_min = rMDI_32b -! Give some initial rMDI values in case they are never set - distance(:,:)=rMDI_32b found=.FALSE. @@ -666,58 +658,42 @@ FUNCTION f_shum_spiral_arg32 & east=-1_int32 ! Find how many points can go south and be inside search_dist - DO j = 1, points_phi - IF (unres_j-j > 0) THEN - tempdist=calc_distance(lats(unres_j),lons(unres_i), & - lats(unres_j-j),lons(unres_i), & - planet_radius) - IF (tempdist > search_dist) THEN - south=j-1_int32 - EXIT - END IF - ELSE + DO j = 1, unres_j-1_int32 + tempdist=calc_distance(lats(unres_j),lons(unres_i), & + lats(unres_j-j),lons(unres_i), & + planet_radius) + IF (tempdist > search_dist) THEN + south=j-1_int32 EXIT END IF END DO ! Find how many points can go north and be inside search_dist - DO j = 1, points_phi - IF (unres_j+j <= points_phi) THEN - tempdist=calc_distance(lats(unres_j),lons(unres_i), & - lats(unres_j+j),lons(unres_i), & - planet_radius) - IF (tempdist > search_dist) THEN - north=j-1_int32 - EXIT - END IF - ELSE + DO j = 1, points_phi-unres_j + tempdist=calc_distance(lats(unres_j),lons(unres_i), & + lats(unres_j+j),lons(unres_i), & + planet_radius) + IF (tempdist > search_dist) THEN + north=j-1_int32 EXIT END IF END DO ! Find how many points can go west and be inside search_dist - DO i = 1, points_lambda - IF (unres_i-i > 0) THEN - tempdist=calc_distance(lats(unres_j),lons(unres_i), & - lats(unres_j),lons(unres_i-i), & - planet_radius) - IF (tempdist > search_dist) THEN - west=i-1_int32 - EXIT - END IF - ELSE + DO i = 1, unres_i-1_int32 + tempdist=calc_distance(lats(unres_j),lons(unres_i), & + lats(unres_j),lons(unres_i-i), & + planet_radius) + IF (tempdist > search_dist) THEN + west=i-1_int32 EXIT END IF END DO ! Find how many points can go east and be inside search_dist - DO i = 1, points_lambda - IF (unres_i+i <= points_lambda) THEN - tempdist=calc_distance(lats(unres_j),lons(unres_i), & - lats(unres_j),lons(unres_i+i), & - planet_radius) - IF (tempdist > search_dist) THEN - east=i-1_int32 - EXIT - END IF - ELSE + DO i = 1, points_lambda-unres_i + tempdist=calc_distance(lats(unres_j),lons(unres_i), & + lats(unres_j),lons(unres_i+i), & + planet_radius) + IF (tempdist > search_dist) THEN + east=i-1_int32 EXIT END IF END DO @@ -749,6 +725,13 @@ FUNCTION f_shum_spiral_arg32 & ! Calculate distance from point distance(i,j) = calc_distance(lats(unres_j), lons(unres_i), & lats(j), lons(i), planet_radius) + END IF + END DO + + DO i = 1, points_lambda + l = i+(j - 1_int32)*points_lambda + ! Check to see if it is a resolved point + IF (.NOT. unres_mask(l)) THEN ! Same type of point. IF (lsm(l) .EQV. is_land_field) THEN ! If current distance is less than any previous min store it. @@ -811,6 +794,13 @@ FUNCTION f_shum_spiral_arg32 & ! Calculate distance from point distance(i,j) = calc_distance(lats(unres_j),lons(unres_i), & lats(j),lons(i), planet_radius) + END IF + END DO + + DO i = unres_i-west, unres_i+east + l = i+(j - 1_int32)*points_lambda + ! Check to see if it is a resolved point + IF (.NOT. unres_mask(l)) THEN ! Same type of point and distance less than maximum distance. IF ((lsm(l) .EQV. is_land_field) .AND. & (distance(i,j) < search_dist)) THEN