diff --git a/components/homme/src/share/cxx/utilities/InternalDiagnostics.cpp b/components/homme/src/share/cxx/utilities/InternalDiagnostics.cpp index 3249966c0b5e..ffc9657f4108 100644 --- a/components/homme/src/share/cxx/utilities/InternalDiagnostics.cpp +++ b/components/homme/src/share/cxx/utilities/InternalDiagnostics.cpp @@ -13,6 +13,8 @@ #include "TimeLevel.hpp" #include "mpi/Comm.hpp" +#include + namespace Homme { namespace { @@ -50,10 +52,10 @@ void print_global_state_hash (const std::string& label) { all_reduce_HashType(comm.mpi_comm(), accum, gaccum, 5); if (comm.root()) { for (int i = 0; i < NUM_TIME_LEVELS; ++i) - fprintf(stderr, "hxxhash> %14d %1d %16lx (E %s)\n", + fprintf(stderr, "hxxhash> %14d %1d %16" PRIx64 " (E %s)\n", tl.nstep, i, gaccum[i], label.c_str()); for (int i = 0; i < Q_NUM_TIME_LEVELS; ++i) - fprintf(stderr, "hxxhash> %14d %1d %16lx (T %s)\n", + fprintf(stderr, "hxxhash> %14d %1d %16" PRIx64 " (T %s)\n", tl.nstep, i, gaccum[NUM_TIME_LEVELS+i], label.c_str()); } } diff --git a/components/homme/src/share/gllfvremap_mod.F90 b/components/homme/src/share/gllfvremap_mod.F90 index e0e0fa6c4daa..e8013df217b1 100644 --- a/components/homme/src/share/gllfvremap_mod.F90 +++ b/components/homme/src/share/gllfvremap_mod.F90 @@ -2088,9 +2088,14 @@ subroutine gfr_f_get_latlon(ie, i, j, lat, lon) type (spherical_polar_t) :: p - p = change_coordinates(gfr%center_f(i,j,ie)) - lat = p%lat - lon = p%lon + if (gfr%is_planar) then + lon = gfr%center_f(i,j,ie)%x + lat = gfr%center_f(i,j,ie)%y + else + p = change_coordinates(gfr%center_f(i,j,ie)) + lat = p%lat + lon = p%lon + end if end subroutine gfr_f_get_latlon subroutine gfr_f_get_corner_latlon(ie, i, j, c, lat, lon) @@ -2103,9 +2108,14 @@ subroutine gfr_f_get_corner_latlon(ie, i, j, c, lat, lon) type (spherical_polar_t) :: p - p = change_coordinates(gfr%corners_f(c,i,j,ie)) - lat = p%lat - lon = p%lon + if (gfr%is_planar) then + lon = gfr%corners_f(c,i,j,ie)%x + lat = gfr%corners_f(c,i,j,ie)%y + else + p = change_coordinates(gfr%corners_f(c,i,j,ie)) + lat = p%lat + lon = p%lon + end if end subroutine gfr_f_get_corner_latlon subroutine gfr_f_get_cartesian3d(ie, i, j, p) @@ -2318,6 +2328,7 @@ subroutine set_ps_Q(elem, nets, nete, timeidx, qidx, nlev) ! Make up a test function for use in unit tests. use coordinate_systems_mod, only: cartesian3D_t, change_coordinates + use physical_constants, only: dd_pi, Lx, Ly type (element_t), intent(inout) :: elem(:) integer, intent(in) :: nets, nete, timeidx, qidx, nlev @@ -2329,7 +2340,13 @@ subroutine set_ps_Q(elem, nets, nete, timeidx, qidx, nlev) do ie = nets, nete do j = 1,np do i = 1,np - p = change_coordinates(elem(ie)%spherep(i,j)) + if (gfr%is_planar) then + p%x = elem(ie)%spherep(i,j)%lon*2*dd_pi/Lx + p%y = elem(ie)%spherep(i,j)%lat*2*dd_pi/Ly + p%z = 0 + else + p = change_coordinates(elem(ie)%spherep(i,j)) + end if elem(ie)%state%ps_v(i,j,timeidx) = & 1.0d3*(1 + 0.05*sin(2*p%x+0.5)*sin(p%y+1.5)*sin(3*p%z+2.5)) q = 0.5*(1 + sin(3*p%x)*sin(3*p%y)*sin(4*p%z)) @@ -2609,7 +2626,8 @@ function check(par, dom_mt, gfr, elem, verbose) result(nerr) logical, intent(in) :: verbose real(kind=real_kind) :: a, b, rd, x, y, f0(np*np), f1(np*np), g(np,np), & - wf(np*np), wg(np,np), qmin, qmax, qmin1, qmax1 + wf(np*np), wg(np,np), qmin, qmax, qmin1, qmax1, lat, lon, latext(2), & + lonext(2), tol integer :: nf, nf2, ie, i, j, k, iremap, info, ilimit, it real(kind=real_kind), allocatable :: Qdp_fv(:,:), ps_v_fv(:,:), & qmins(:,:,:), qmaxs(:,:,:) @@ -2669,6 +2687,31 @@ function check(par, dom_mt, gfr, elem, verbose) result(nerr) nerr = nerr + 1 write(iulog,*) 'gfr> Dinv', ie, rd end if + if (gfr%is_planar) then + lonext(1) = minval(elem(ie)%spherep(:,:)%lon) + lonext(2) = maxval(elem(ie)%spherep(:,:)%lon) + latext(1) = minval(elem(ie)%spherep(:,:)%lat) + latext(2) = maxval(elem(ie)%spherep(:,:)%lat) + tol = max(lonext(2) - lonext(1), latext(2) - latext(1))*1.e4*eps + do i = 1, nf + do j = 1, nf + call gfr_f_get_latlon(ie, i, j, lat, lon) + if ( lon < lonext(1) .or. lon > lonext(2) .or. & + lat < latext(1) .or. lat > latext(2)) then + write(iulog,*) 'gfr> planar latlon ctr', ie, latext, lonext, lat, lon + nerr = nerr + 1 + end if + do k = 1, 4 + call gfr_f_get_corner_latlon(ie, i, j, k, lat, lon) + if ( lon < lonext(1) - tol .or. lon > lonext(2) + tol .or. & + lat < latext(1) - tol .or. lat > latext(2) + tol) then + write(iulog,*) 'gfr> planar latlon crn', ie, latext, lonext, lat, lon + nerr = nerr + 1 + end if + end do + end do + end do + end if ! Check that FV -> GLL -> FV recovers the original FV values exactly ! (with no DSS and no limiter). @@ -2839,8 +2882,8 @@ function test_sphere2ref() result(nerr) use kinds, only: iulog type (cartesian3D_t) :: corners(4), sphere - real (real_kind) :: refin(2), refout(2), err - integer :: i, j, n, nerr + real (real_kind) :: refin(2), refout(2), err, tol + integer :: i, j, n, trial, nerr nerr = 0 @@ -2854,18 +2897,28 @@ function test_sphere2ref() result(nerr) do i = 1,n refin(1) = -1 + (1.0_real_kind/(n-1))*(i-1) do j = 1,n - refin(2) = -1 + (1.0_real_kind/(n-1))*(j-1) - call ref2sphere(corners, refin(1), refin(2), sphere) - call sphere2ref(corners, sphere, refout(1), refout(2)) - err = abs(refin(1) - refout(1)) + abs(refin(2) - refout(2)) - if (err > 15*eps .or. & - maxval(abs(refout)) > 1 + 5*eps .or. & - any(refout /= refout)) then - write(iulog,*) refin(1), refin(2) - write(iulog,*) refout(1), refout(2) - write(iulog,*) err - nerr = nerr + 1 - end if + do trial = 1, 2 + refin(2) = -1 + (1.0_real_kind/(n-1))*(j-1) + call ref2sphere(corners, refin(1), refin(2), sphere) + if (trial == 2) then + sphere%x = 0.9d0*sphere%x + sphere%y = 0.9d0*sphere%y + sphere%z = 0.9d0*sphere%z + end if + call sphere2ref(corners, sphere, refout(1), refout(2)) + err = abs(refin(1) - refout(1)) + abs(refin(2) - refout(2)) + tol = 15*eps + if (trial == 2) tol = 1e-6 + if (err > tol .or. & + maxval(abs(refout)) > 1 + 5*eps .or. & + any(refout /= refout)) then + write(iulog,*) 'gfr> test_sphere2ref trial', trial + write(iulog,*) refin(1), refin(2) + write(iulog,*) refout(1), refout(2) + write(iulog,*) err + nerr = nerr + 1 + end if + end do end do end do if (nerr /= 0) write(iulog,*) 'test_sphere2ref FAILED'