From 23bbfcd6013447e248be5e64e834b5b98dc075cf Mon Sep 17 00:00:00 2001 From: Bin Liu Date: Tue, 17 Jan 2023 15:56:59 -0500 Subject: [PATCH] Update submodule UPP and some HAFS moving-nesting related fixes (#613) * Hardcoded coarse terrain for moving nest with terrain_smoother=4 * Update HAFSv1 CCPP suites to use the unified ugwp scheme. * Update submodule upp. * Removed debug print statements for terrain_smoother * Fixed bug in calculation of nest halo weights identified by Biju Thomas in DDEBUG=ON builds. * Updated fix to moving nest weight calculations. * Use ind arrays for nest weight calculations. Ensures correct indexing. * Clean up for the fix of using ind arrays for moving-nest weight calculations (from @wramstrom). Co-authored-by: William Ramstrom --- .../suites/suite_FV3_HAFS_v1_gfdlmp_tedmf.xml | 4 +- .../suite_FV3_HAFS_v1_gfdlmp_tedmf_nonsst.xml | 4 +- ccpp/suites/suite_FV3_HAFS_v1_thompson.xml | 4 +- .../suite_FV3_HAFS_v1_thompson_noahmp.xml | 4 +- ...ite_FV3_HAFS_v1_thompson_noahmp_nonsst.xml | 4 +- .../suite_FV3_HAFS_v1_thompson_nonsst.xml | 4 +- ...uite_FV3_HAFS_v1_thompson_tedmf_gfdlsf.xml | 4 +- moving_nest/fv_moving_nest.F90 | 303 +++++++++++++++--- moving_nest/fv_moving_nest_main.F90 | 32 +- moving_nest/fv_moving_nest_types.F90 | 4 +- moving_nest/fv_moving_nest_utils.F90 | 2 +- upp | 2 +- 12 files changed, 307 insertions(+), 64 deletions(-) diff --git a/ccpp/suites/suite_FV3_HAFS_v1_gfdlmp_tedmf.xml b/ccpp/suites/suite_FV3_HAFS_v1_gfdlmp_tedmf.xml index 1760afe97..3dc3c8d54 100644 --- a/ccpp/suites/suite_FV3_HAFS_v1_gfdlmp_tedmf.xml +++ b/ccpp/suites/suite_FV3_HAFS_v1_gfdlmp_tedmf.xml @@ -62,8 +62,8 @@ satmedmfvdifq GFS_PBL_generic_post GFS_GWD_generic_pre - cires_ugwp - cires_ugwp_post + unified_ugwp + unified_ugwp_post GFS_GWD_generic_post GFS_suite_stateout_update ozphys_2015 diff --git a/ccpp/suites/suite_FV3_HAFS_v1_gfdlmp_tedmf_nonsst.xml b/ccpp/suites/suite_FV3_HAFS_v1_gfdlmp_tedmf_nonsst.xml index 90693645c..254df77e0 100644 --- a/ccpp/suites/suite_FV3_HAFS_v1_gfdlmp_tedmf_nonsst.xml +++ b/ccpp/suites/suite_FV3_HAFS_v1_gfdlmp_tedmf_nonsst.xml @@ -60,8 +60,8 @@ satmedmfvdifq GFS_PBL_generic_post GFS_GWD_generic_pre - cires_ugwp - cires_ugwp_post + unified_ugwp + unified_ugwp_post GFS_GWD_generic_post GFS_suite_stateout_update ozphys_2015 diff --git a/ccpp/suites/suite_FV3_HAFS_v1_thompson.xml b/ccpp/suites/suite_FV3_HAFS_v1_thompson.xml index 0e7127bb4..746c8ceb7 100644 --- a/ccpp/suites/suite_FV3_HAFS_v1_thompson.xml +++ b/ccpp/suites/suite_FV3_HAFS_v1_thompson.xml @@ -57,8 +57,8 @@ satmedmfvdifq GFS_PBL_generic_post GFS_GWD_generic_pre - cires_ugwp - cires_ugwp_post + unified_ugwp + unified_ugwp_post GFS_GWD_generic_post GFS_suite_stateout_update ozphys_2015 diff --git a/ccpp/suites/suite_FV3_HAFS_v1_thompson_noahmp.xml b/ccpp/suites/suite_FV3_HAFS_v1_thompson_noahmp.xml index 8d3bab273..2022ef4d6 100644 --- a/ccpp/suites/suite_FV3_HAFS_v1_thompson_noahmp.xml +++ b/ccpp/suites/suite_FV3_HAFS_v1_thompson_noahmp.xml @@ -57,8 +57,8 @@ satmedmfvdifq GFS_PBL_generic_post GFS_GWD_generic_pre - cires_ugwp - cires_ugwp_post + unified_ugwp + unified_ugwp_post GFS_GWD_generic_post GFS_suite_stateout_update ozphys_2015 diff --git a/ccpp/suites/suite_FV3_HAFS_v1_thompson_noahmp_nonsst.xml b/ccpp/suites/suite_FV3_HAFS_v1_thompson_noahmp_nonsst.xml index b61422021..f66543c80 100644 --- a/ccpp/suites/suite_FV3_HAFS_v1_thompson_noahmp_nonsst.xml +++ b/ccpp/suites/suite_FV3_HAFS_v1_thompson_noahmp_nonsst.xml @@ -55,8 +55,8 @@ satmedmfvdifq GFS_PBL_generic_post GFS_GWD_generic_pre - cires_ugwp - cires_ugwp_post + unified_ugwp + unified_ugwp_post GFS_GWD_generic_post GFS_suite_stateout_update ozphys_2015 diff --git a/ccpp/suites/suite_FV3_HAFS_v1_thompson_nonsst.xml b/ccpp/suites/suite_FV3_HAFS_v1_thompson_nonsst.xml index 0cb393317..665a0b7e2 100644 --- a/ccpp/suites/suite_FV3_HAFS_v1_thompson_nonsst.xml +++ b/ccpp/suites/suite_FV3_HAFS_v1_thompson_nonsst.xml @@ -55,8 +55,8 @@ satmedmfvdifq GFS_PBL_generic_post GFS_GWD_generic_pre - cires_ugwp - cires_ugwp_post + unified_ugwp + unified_ugwp_post GFS_GWD_generic_post GFS_suite_stateout_update ozphys_2015 diff --git a/ccpp/suites/suite_FV3_HAFS_v1_thompson_tedmf_gfdlsf.xml b/ccpp/suites/suite_FV3_HAFS_v1_thompson_tedmf_gfdlsf.xml index e68c9e687..a5db1110b 100644 --- a/ccpp/suites/suite_FV3_HAFS_v1_thompson_tedmf_gfdlsf.xml +++ b/ccpp/suites/suite_FV3_HAFS_v1_thompson_tedmf_gfdlsf.xml @@ -57,8 +57,8 @@ satmedmfvdifq GFS_PBL_generic_post GFS_GWD_generic_pre - cires_ugwp - cires_ugwp_post + unified_ugwp + unified_ugwp_post GFS_GWD_generic_post GFS_suite_stateout_update ozphys_2015 diff --git a/moving_nest/fv_moving_nest.F90 b/moving_nest/fv_moving_nest.F90 index 4e2ab8b59..009f04b87 100644 --- a/moving_nest/fv_moving_nest.F90 +++ b/moving_nest/fv_moving_nest.F90 @@ -54,7 +54,7 @@ module fv_moving_nest_mod use block_control_mod, only : block_control_type use fms_mod, only : mpp_clock_id, mpp_clock_begin, mpp_clock_end, CLOCK_ROUTINE, clock_flag_default, CLOCK_SUBCOMPONENT - use mpp_mod, only : mpp_pe, mpp_sync, mpp_sync_self, mpp_send, mpp_error, NOTE, FATAL + use mpp_mod, only : mpp_pe, mpp_sync, mpp_sync_self, mpp_send, mpp_error, NOTE, FATAL, WARNING use mpp_domains_mod, only : mpp_update_domains, mpp_get_data_domain, mpp_get_global_domain use mpp_domains_mod, only : mpp_define_nest_domains, mpp_shift_nest_domains, nest_domain_type, domain2d use mpp_domains_mod, only : mpp_get_C2F_index, mpp_update_nest_fine @@ -261,6 +261,12 @@ subroutine mn_prog_fill_nest_halos_from_parent(Atm, n, child_grid_num, is_fine_p call fill_nest_halos_from_parent("q", Atm(n)%q, interp_type, Atm(child_grid_num)%neststruct%wt_h, & Atm(child_grid_num)%neststruct%ind_h, x_refine, y_refine, is_fine_pe, nest_domain, position, nz) + ! Interpolate terrain from coarse grid + if (Moving_nest(n)%mn_flag%terrain_smoother .eq. 4) then + call fill_nest_halos_from_parent("phis", Atm(n)%phis, interp_type, Atm(child_grid_num)%neststruct%wt_h, & + Atm(child_grid_num)%neststruct%ind_h, x_refine, y_refine, is_fine_pe, nest_domain, position) + endif + ! Move the A-grid winds. TODO consider recomputing them from D grid instead call fill_nest_halos_from_parent("ua", Atm(n)%ua, interp_type, Atm(child_grid_num)%neststruct%wt_h, & Atm(child_grid_num)%neststruct%ind_h, x_refine, y_refine, is_fine_pe, nest_domain, position, nz) @@ -350,6 +356,8 @@ subroutine mn_prog_fill_intern_nest_halos(Atm, domain_fine, is_fine_pe) integer :: this_pe type(fv_moving_nest_prog_type), pointer :: mn_prog + integer :: child_grid_num = 2 + mn_prog => Moving_nest(2)%mn_prog ! TODO allow nest number to vary this_pe = mpp_pe() @@ -360,6 +368,10 @@ subroutine mn_prog_fill_intern_nest_halos(Atm, domain_fine, is_fine_pe) call mn_var_fill_intern_nest_halos(Atm%delp, domain_fine, is_fine_pe) call mn_var_fill_intern_nest_halos(mn_prog%delz, domain_fine, is_fine_pe) + if (Moving_nest(child_grid_num)%mn_flag%terrain_smoother .eq. 4) then + call mn_var_fill_intern_nest_halos(Atm%phis, domain_fine, is_fine_pe) + endif + call mn_var_fill_intern_nest_halos(Atm%ua, domain_fine, is_fine_pe) call mn_var_fill_intern_nest_halos(Atm%va, domain_fine, is_fine_pe) @@ -596,14 +608,27 @@ subroutine mn_latlon_load_parent(surface_dir, Atm, n, parent_tile, delta_i_c, de if (first_nest_move) then if (use_timers) call mpp_clock_begin (id_load1) + parent_geo%nxp = Atm(1)%npx + parent_geo%nyp = Atm(1)%npy + + parent_geo%nx = parent_geo%nxp - 1 + parent_geo%ny = parent_geo%nyp - 1 + call mn_static_filename(surface_dir, parent_tile, 'grid', 1, grid_filename) - call load_nest_latlons_from_nc(grid_filename, Atm(1)%npx, Atm(1)%npy, 1, pelist, & + call load_nest_latlons_from_nc(grid_filename, parent_geo%nxp, parent_geo%nyp, 1, pelist, & parent_geo, p_istart_fine, p_iend_fine, p_jstart_fine, p_jend_fine) + ! These are saved between timesteps in fv_moving_nest_main.F90 - allocate(p_grid(1:parent_geo%nxp, 1:parent_geo%nyp,2)) - allocate(p_grid_u(1:parent_geo%nxp, 1:parent_geo%nyp+1,2)) - allocate(p_grid_v(1:parent_geo%nxp+1, 1:parent_geo%nyp,2)) + allocate(p_grid(parent_geo%nx, parent_geo%ny,2)) + allocate(p_grid_u(parent_geo%nx, parent_geo%ny+1,2)) + allocate(p_grid_v(parent_geo%nx+1, parent_geo%ny,2)) + + p_grid = 0.0 + p_grid_u = 0.0 + p_grid_v = 0.0 + + ! These are big (parent grid size), and do not change during the model integration. call assign_p_grids(parent_geo, p_grid, position) @@ -616,11 +641,6 @@ subroutine mn_latlon_load_parent(surface_dir, Atm, n, parent_tile, delta_i_c, de if (use_timers) call mpp_clock_begin (id_load2) - parent_geo%nxp = Atm(1)%npx - parent_geo%nyp = Atm(1)%npy - - parent_geo%nx = Atm(1)%npx - 1 - parent_geo%ny = Atm(1)%npy - 1 !=========================================================== ! Begin tile_geo per PE. @@ -772,6 +792,8 @@ subroutine mn_latlon_read_hires_parent(npx, npy, refine, pelist, fp_super_tile_g integer :: fp_super_istart_fine, fp_super_jstart_fine,fp_super_iend_fine, fp_super_jend_fine character(len=256) :: grid_filename + integer :: i, j + call mn_static_filename(surface_dir, parent_tile, 'grid', refine, grid_filename) call load_nest_latlons_from_nc(trim(grid_filename), npx, npy, refine, pelist, & @@ -910,7 +932,7 @@ end subroutine mn_static_read_hires_r8 !>@brief The subroutine 'mn_meta_recalc' recalculates nest halo weights subroutine mn_meta_recalc( delta_i_c, delta_j_c, x_refine, y_refine, tile_geo, parent_geo, fp_super_tile_geo, & - is_fine_pe, nest_domain, position, p_grid, n_grid, wt, istart_coarse, jstart_coarse) + is_fine_pe, nest_domain, position, p_grid, n_grid, wt, istart_coarse, jstart_coarse, ind) integer, intent(in) :: delta_i_c, delta_j_c !< Nest motion in delta i,j integer, intent(in) :: x_refine, y_refine !< Nest refinement type(grid_geometry), intent(inout) :: tile_geo, parent_geo, fp_super_tile_geo !< tile geometries @@ -921,8 +943,10 @@ subroutine mn_meta_recalc( delta_i_c, delta_j_c, x_refine, y_refine, tile_geo, p real, allocatable, intent(inout) :: wt(:,:,:) !< Interpolation weights integer, intent(inout) :: position !< Stagger integer, intent(in) :: istart_coarse, jstart_coarse !< Initian nest offsets + integer, allocatable, intent(in) :: ind(:,:,:) type(bbox) :: wt_fine, wt_coarse + integer :: istag, jstag integer :: this_pe this_pe = mpp_pe() @@ -936,17 +960,41 @@ subroutine mn_meta_recalc( delta_i_c, delta_j_c, x_refine, y_refine, tile_geo, p !! !!=========================================================== + if (position .eq. CENTER) then + istag = 0 + jstag = 0 + elseif (position .eq. NORTH) then + ! for p_grid_u + istag = 1 + jstag = 0 + + !! This aligns with boundary.F90 settings + !!istag = 0 + !!jstag = 1 + + elseif (position .eq. EAST) then + ! for p_grid_v + istag = 0 + jstag = 1 + + !! This aligns with boundary.F90 settings + !istag = 1 + !jstag = 0 + + endif + + call bbox_get_C2F_index(nest_domain, wt_fine, wt_coarse, EAST, position) - call calc_nest_halo_weights(wt_fine, wt_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine) + call calc_nest_halo_weights(wt_fine, wt_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine, istag, jstag, ind) call bbox_get_C2F_index(nest_domain, wt_fine, wt_coarse, WEST, position) - call calc_nest_halo_weights(wt_fine, wt_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine) + call calc_nest_halo_weights(wt_fine, wt_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine, istag, jstag, ind) call bbox_get_C2F_index(nest_domain, wt_fine, wt_coarse, NORTH, position) - call calc_nest_halo_weights(wt_fine, wt_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine) + call calc_nest_halo_weights(wt_fine, wt_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine, istag, jstag, ind) call bbox_get_C2F_index(nest_domain, wt_fine, wt_coarse, SOUTH, position) - call calc_nest_halo_weights(wt_fine, wt_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine) + call calc_nest_halo_weights(wt_fine, wt_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine, istag, jstag, ind) endif @@ -1034,6 +1082,11 @@ subroutine mn_prog_shift_data(Atm, n, child_grid_num, wt_h, wt_u, wt_v, & call mn_var_shift_data(mn_prog%delz, interp_type, wt_h, Atm(child_grid_num)%neststruct%ind_h, & delta_i_c, delta_j_c, x_refine, y_refine, is_fine_pe, nest_domain, position, nz) + if (Moving_nest(n)%mn_flag%terrain_smoother .eq. 4) then + call mn_var_shift_data(Atm(n)%phis, interp_type, wt_h, Atm(child_grid_num)%neststruct%ind_h, & + delta_i_c, delta_j_c, x_refine, y_refine, is_fine_pe, nest_domain, position) + endif + call mn_var_shift_data(Atm(n)%ua, interp_type, wt_h, Atm(child_grid_num)%neststruct%ind_h, & delta_i_c, delta_j_c, x_refine, y_refine, is_fine_pe, nest_domain, position, nz) @@ -2371,36 +2424,123 @@ subroutine assign_p_grids(parent_geo, p_grid, position) integer :: i,j + integer :: num_zeros, num_vals, full_zeros + + num_zeros = 0 + full_zeros = 0 + num_vals = 0 + if (position == CENTER) then - do j = 1, parent_geo%ny - do i = 1, parent_geo%nx + do j = 1, ubound(p_grid,2) + do i = 1, ubound(p_grid,1) ! centered grid version - p_grid(i, j, 1) = parent_geo%lons(2*i, 2*j) - p_grid(i, j, 2) = parent_geo%lats(2*i, 2*j) + p_grid(i, j, :) = 0.0 + + if (2*i .gt. ubound(parent_geo%lons,1)) then + print '("[ERROR] WDR PG_CLONi npe=",I0," 2*i=",I0," ubound=",I0)', mpp_pe(), 2*i, ubound(parent_geo%lons,1) + elseif (2*j .gt. ubound(parent_geo%lons,2)) then + print '("[ERROR] WDR PG_CLONj npe=",I0," 2*j=",I0," ubound=",I0)', mpp_pe(), 2*j, ubound(parent_geo%lons,2) + else + p_grid(i, j, 1) = parent_geo%lons(2*i, 2*j) + p_grid(i, j, 2) = parent_geo%lats(2*i, 2*j) + endif enddo enddo + + + do i = 1, ubound(p_grid,1) + do j = 1, ubound(p_grid,2) + + if (p_grid(i,j,1) .eq. 0.0) then + num_zeros = num_zeros + 1 + if (p_grid(i,j,2) .eq. 0.0) then + full_zeros = full_zeros + 1 + print '("[INFO] WDR set p_grid FULL_ZERO npe=",I0," i=",I0," j=",I0)', mpp_pe(), i, j + endif + endif + if (p_grid(i,j,1) .ne. 0.0) num_vals = num_vals + 1 + enddo + enddo + + if (num_zeros .gt. 0) print '("[INFO] WDR set p_grid npe=",I0," num_zeros=",I0," full_zeros=",I0," num_vals=",I0" nxp=",I0," nyp=",I0," parent_geo%lats(",I0,",",I0,")"," p_grid(",I0,",",I0,",2)")', mpp_pe(), num_zeros, full_zeros, num_vals, parent_geo%nxp, parent_geo%nyp, ubound(parent_geo%lats,1), ubound(parent_geo%lats,2), ubound(p_grid,1), ubound(p_grid,2) + ! u(npx, npy+1) elseif (position == NORTH) then ! u wind on D-stagger - do j = 1, parent_geo%ny - do i = 1, parent_geo%nx - ! centered grid version - p_grid(i, j, 1) = parent_geo%lons(2*i, 2*j-1) + do j = 1, ubound(p_grid,2) + do i = 1, ubound(p_grid,1) + ! u grid version + p_grid(i, j, :) = 0.0 + if (2*i .gt. ubound(parent_geo%lons,1)) then + print '("[ERROR] WDR PG_ULONi npe=",I0," 2*i=",I0," ubound=",I0)', mpp_pe(), 2*i, ubound(parent_geo%lons,1) + elseif (2*j-1 .gt. ubound(parent_geo%lons,2)) then + print '("[ERROR] WDR PG_ULONj npe=",I0," 2*j-1=",I0," ubound=",I0)', mpp_pe(), 2*j-1, ubound(parent_geo%lons,2) + else + + ! This seems correct + p_grid(i, j, 1) = parent_geo%lons(2*i, 2*j-1) p_grid(i, j, 2) = parent_geo%lats(2*i, 2*j-1) - enddo + endif enddo + enddo + + + do i = 1, ubound(p_grid,1) + do j = 1, ubound(p_grid,2) + + if (p_grid(i,j,1) .eq. 0.0) then + num_zeros = num_zeros + 1 + if (p_grid(i,j,2) .eq. 0.0) then + full_zeros = full_zeros + 1 + print '("[INFO] WDR set p_grid_u FULL_ZERO npe=",I0," i=",I0," j=",I0)', mpp_pe(), i, j + endif + endif + if (p_grid(i,j,1) .ne. 0.0) num_vals = num_vals + 1 + enddo + enddo + + if (num_zeros .gt. 0) print '("[INFO] WDR set p_grid_u npe=",I0," num_zeros=",I0," full_zeros=",I0," num_vals=",I0" nxp=",I0," nyp=",I0," parent_geo%lats(",I0,",",I0,")"," p_grid(",I0,",",I0,",2)")', mpp_pe(), num_zeros, full_zeros, num_vals, parent_geo%nxp, parent_geo%nyp, ubound(parent_geo%lats,1), ubound(parent_geo%lats,2), ubound(p_grid,1), ubound(p_grid,2) + ! v(npx+1, npy) elseif (position == EAST) then ! v wind on D-stagger - do j = 1, parent_geo%ny - do i = 1, parent_geo%nx - ! centered grid version - p_grid(i, j, 1) = parent_geo%lons(2*i-1, 2*j) - p_grid(i, j, 2) = parent_geo%lats(2*i-1, 2*j) + do j = 1, ubound(p_grid,2) + do i = 1, ubound(p_grid,1) + ! v grid version + p_grid(i, j, :) = 0.0 + if (2*i-1 .gt. ubound(parent_geo%lons,1)) then + print '("[ERROR] WDR PG_VLONi npe=",I0," 2*i-1=",I0," ubound=",I0)', mpp_pe(), 2*i-1, ubound(parent_geo%lons,1) + elseif (2*j .gt. ubound(parent_geo%lons,2)) then + print '("[ERROR] WDR PG_VLONj npe=",I0," 2*j=",I0," ubound=",I0)', mpp_pe(), 2*j, ubound(parent_geo%lons,2) + else + ! This seems correct + p_grid(i, j, 1) = parent_geo%lons(2*i-1, 2*j) + p_grid(i, j, 2) = parent_geo%lats(2*i-1, 2*j) + endif + enddo + enddo + + + do i = 1, ubound(p_grid,1) + do j = 1, ubound(p_grid,2) + + if (p_grid(i,j,1) .eq. 0.0) then + num_zeros = num_zeros + 1 + if (p_grid(i,j,2) .eq. 0.0) then + full_zeros = full_zeros + 1 + print '("[INFO] WDR set p_grid_v FULL_ZERO npe=",I0," i=",I0," j=",I0)', mpp_pe(), i, j + endif + endif + if (p_grid(i,j,1) .ne. 0.0) num_vals = num_vals + 1 enddo enddo + + if (num_zeros .gt. 0) print '("[INFO] WDR set p_grid_v npe=",I0," num_zeros=",I0," full_zeros=",I0," num_vals=",I0" nxp=",I0," nyp=",I0," parent_geo%lats(",I0,",",I0,")"," p_grid(",I0,",",I0,",2)")', mpp_pe(), num_zeros, full_zeros, num_vals, parent_geo%nxp, parent_geo%nyp, ubound(parent_geo%lats,1), ubound(parent_geo%lats,2), ubound(p_grid,1), ubound(p_grid,2) + + + endif - + end subroutine assign_p_grids @@ -2447,28 +2587,65 @@ subroutine assign_n_grids(tile_geo, n_grid, position) end subroutine assign_n_grids + subroutine calc_inside(p_grid, ic, jc, n_grid1, n_grid2, istag, jstag, is_inside, verbose) + real(kind=R_GRID), allocatable, intent(in) :: p_grid(:,:,:) + real(kind=R_GRID), intent(in) :: n_grid1, n_grid2 + integer, intent(in) :: ic, jc, istag, jstag + logical, intent(out) :: is_inside + logical, intent(in) :: verbose + + + real(kind=R_GRID) :: max1, max2, min1, min2, eps + + max1 = max(p_grid(ic,jc,1), p_grid(ic,jc+1,1), p_grid(ic,jc+1,1), p_grid(ic+1,jc+1,1), p_grid(ic+1,jc+1,1), p_grid(ic+1,jc,1)) + max2 = max(p_grid(ic,jc,2), p_grid(ic,jc+1,2), p_grid(ic,jc+1,2), p_grid(ic+1,jc+1,2), p_grid(ic+1,jc+1,2), p_grid(ic+1,jc,2)) + + min1 = min(p_grid(ic,jc,1), p_grid(ic,jc+1,1), p_grid(ic,jc+1,1), p_grid(ic+1,jc+1,1), p_grid(ic+1,jc+1,1), p_grid(ic+1,jc,1)) + min2 = min(p_grid(ic,jc,2), p_grid(ic,jc+1,2), p_grid(ic,jc+1,2), p_grid(ic+1,jc+1,2), p_grid(ic+1,jc+1,2), p_grid(ic+1,jc,2)) + + is_inside = .False. + + eps = 0.00001 + !eps = 0.000001 + + if (n_grid1 .le. max1+eps .and. n_grid1 .ge. min1-eps) then + if (n_grid2 .le. max2+eps .and. n_grid2 .ge. min2-eps) then + is_inside = .True. + !if (verbose) print '("[INFO] WDR is_inside TRUE npe=",I0," ic=",I0," jc=",I0," n_grid1=",F12.8," min1=",F12.8," max1=",F12.8," n_grid2=",F12.8," min2=",F12.8," max2=", F12.8," p_grid(",I0,",",I0,",2) istag=",I0," jstag=",I0)', mpp_pe(), ic, jc, n_grid1, min1, max1, n_grid2, min2, max2, ubound(p_grid,1), ubound(p_grid,2), istag, jstag + endif + endif + + if (verbose .and. .not. is_inside) then + print '("[INFO] WDR is_inside FALSE npe=",I0," ic=",I0," jc=",I0," n_grid1=",F12.8," min1=",F12.8," max1=",F12.8," n_grid2=",F12.8," min2=",F12.8," max2=", F10.4," p_grid(",I0,",",I0,",2) istag=",I0," jstag=",I0)', mpp_pe(), ic, jc, n_grid1, min1, max1, n_grid2, min2, max2, ubound(p_grid,1), ubound(p_grid,2), istag, jstag + endif + end subroutine calc_inside + !>@brief The subroutine 'calc_nest_halo_weights' calculates the interpolation weights !>@details Computationally demanding; target for optimization after nest moves - subroutine calc_nest_halo_weights(bbox_fine, bbox_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine) + subroutine calc_nest_halo_weights(bbox_fine, bbox_coarse, p_grid, n_grid, wt, istart_coarse, jstart_coarse, x_refine, y_refine, istag, jstag, ind) implicit none type(bbox), intent(in) :: bbox_coarse, bbox_fine !< Bounding boxes of parent and nest real(kind=R_GRID), allocatable, intent(in) :: p_grid(:,:,:), n_grid(:,:,:) !< Latlon rids of parent and nest in radians real, allocatable, intent(inout) :: wt(:,:,:) !< Interpolation weight array integer, intent(in) :: istart_coarse, jstart_coarse, x_refine, y_refine !< Offsets and nest refinements + integer, intent(in) :: istag, jstag !< Staggers + integer, allocatable, intent(in) :: ind(:,:,:) - integer :: i,j, ic, jc + integer :: i, j, k, ic, jc real :: dist1, dist2, dist3, dist4, sum logical :: verbose = .false. !logical :: verbose = .true. - + logical :: is_inside, adjusted integer :: this_pe real(kind=R_GRID) :: pi = 4 * atan(1.0d0) real :: pi180 real :: rad2deg, deg2rad + real :: old_weight(4), diff_weight(4) + integer :: di, dj pi180 = pi / 180.0 deg2rad = pi / 180.0 @@ -2486,25 +2663,75 @@ subroutine calc_nest_halo_weights(bbox_fine, bbox_coarse, p_grid, n_grid, wt, is ! do j = bbox_fine%js, bbox_fine%je - ! F90 integer division truncates - jc = jstart_coarse + (j + y_refine/2 + 1) / y_refine do i = bbox_fine%is, bbox_fine%ie - ic = istart_coarse + (i + x_refine/2 + 1) / x_refine - + jc = ind(i,j,2) ! reset this if the UPDATE code altered it + ic = ind(i,j,1) + if (ic+1 .gt. ubound(p_grid, 1)) print '("[ERROR] WDR CALCWT off end of p_grid i npe=",I0," ic+1=",I0," bound=",I0)', mpp_pe(), ic+1, ubound(p_grid,1) + if (jc+1 .gt. ubound(p_grid, 2)) print '("[ERROR] WDR CALCWT off end of p_grid j npe=",I0," jc+1=",I0," bound=",I0)', mpp_pe(), jc+1, ubound(p_grid,2) + ! dist2side_latlon takes points in longitude-latitude coordinates. dist1 = dist2side_latlon(p_grid(ic,jc,:), p_grid(ic,jc+1,:), n_grid(i,j,:)) dist2 = dist2side_latlon(p_grid(ic,jc+1,:), p_grid(ic+1,jc+1,:), n_grid(i,j,:)) dist3 = dist2side_latlon(p_grid(ic+1,jc+1,:), p_grid(ic+1,jc,:), n_grid(i,j,:)) dist4 = dist2side_latlon(p_grid(ic,jc,:), p_grid(ic+1,jc,:), n_grid(i,j,:)) + call calc_inside(p_grid, ic, jc, n_grid(i,j,1), n_grid(i,j,2), istag, jstag, is_inside, .True.) + +! if (.not. is_inside) then +! adjusted = .False. +! +! do di = -2,2 +! do dj = -2,1 +! if (.not. adjusted) then +! call calc_inside(p_grid, ic+di, jc+dj, n_grid(i,j,1), n_grid(i,j,2), istag, jstag, is_inside, .False.) +! if (is_inside) then +! ic = ic + di +! jc = jc + dj +! +! print '("[INFO] WDR is_inside UPDATED npe=",I0," ic=",I0," jc=",I0," istart_coarse=",I0," jstart_coarse=",I0," i=",I0," j=",I0," di=",I0," dj=",I0," n_grid1=",F12.8," n_grid2=",F12.8," istag=",I0," jstag=",I0)', mpp_pe(), ic, jc, istart_coarse, jstart_coarse, i, j, di, dj, n_grid(i,j,1), n_grid(i,j,2), istag, jstag + +! dist1 = dist2side_latlon(p_grid(ic,jc,:), p_grid(ic,jc+1,:), n_grid(i,j,:)) +! dist2 = dist2side_latlon(p_grid(ic,jc+1,:), p_grid(ic+1,jc+1,:), n_grid(i,j,:)) +! dist3 = dist2side_latlon(p_grid(ic+1,jc+1,:), p_grid(ic+1,jc,:), n_grid(i,j,:)) +! dist4 = dist2side_latlon(p_grid(ic,jc,:), p_grid(ic+1,jc,:), n_grid(i,j,:)) +! +! adjusted = .True. +! endif +! endif +! enddo +! enddo +! if (.not. adjusted) print '("[ERROR] WDR is_inside UPDATE FAILED npe=",I0," i=",I0," j=",I0," ic=",I0," jc=",I0," n_grid1=",F12.8," n_grid2=",F12.8," istag=",I0," jstag=",I0)', mpp_pe(), i, j, ic, jc, n_grid(i,j,1), n_grid(i,j,2), istag, jstag +! +! endif + + old_weight = wt(i,j,:) + wt(i,j,1)=dist2*dist3 ! ic, jc weight wt(i,j,2)=dist3*dist4 ! ic, jc+1 weight wt(i,j,3)=dist4*dist1 ! ic+1, jc+1 weight wt(i,j,4)=dist1*dist2 ! ic+1, jc weight + ! If sum is zero, this is a problem + sum=wt(i,j,1)+wt(i,j,2)+wt(i,j,3)+wt(i,j,4) - wt(i,j,:)=wt(i,j,:)/sum + if (sum .eq. 0.0) then + + call mpp_error(WARNING, "WARNING: calc_nest_halo_weights sum of weights is zero.") + wt(i,j,:) = 0.25 + + else + wt(i,j,:)=wt(i,j,:)/sum + endif + + diff_weight = old_weight - wt(i,j,:) + do k=1,4 + if (abs(diff_weight(k)) .ge. 0.01) then + print '("[WARN] WDR DIFFWT npe=",I0," old_wt=",F10.6," wt(",I0,",",I0,",",I0,")=",F10.6," diff=",F10.6," istag=",I0," jstag=",I0)', & + mpp_pe(), old_weight(k), i, j, k, wt(i,j,k), diff_weight(k), istag, jstag + + endif + enddo enddo enddo endif diff --git a/moving_nest/fv_moving_nest_main.F90 b/moving_nest/fv_moving_nest_main.F90 index c848e313f..4f1ce7aeb 100644 --- a/moving_nest/fv_moving_nest_main.F90 +++ b/moving_nest/fv_moving_nest_main.F90 @@ -682,6 +682,12 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, allocate(wt_v(Atm(child_grid_num)%bd%isd:Atm(child_grid_num)%bd%ied+1, Atm(child_grid_num)%bd%jsd:Atm(child_grid_num)%bd%jed, 4)) wt_v = real_snan + + ! Fill in the local weights with the ones from Atm just to be safe + call fill_weight_grid(wt_h, Atm(n)%neststruct%wt_h) + call fill_weight_grid(wt_u, Atm(n)%neststruct%wt_u) + call fill_weight_grid(wt_v, Atm(n)%neststruct%wt_v) + else allocate(wt_h(1,1,4)) wt_h = 0.0 @@ -912,20 +918,28 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, call mn_reset_phys_latlon(Atm, n, tile_geo, fp_super_tile_geo, Atm_block, IPD_control, IPD_data) if (use_timers) call mpp_clock_end (id_movnest5_2) - if (use_timers) call mpp_clock_begin (id_movnest5_3) + endif !!============================================================================ !! Step 5.2 -- Fill the wt* variables for each stagger !!============================================================================ + call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_h) + call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_u) + call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_v) + call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_b) + + if (is_fine_pe) then + if (use_timers) call mpp_clock_begin (id_movnest5_3) + call mn_meta_recalc( delta_i_c, delta_j_c, x_refine, y_refine, tile_geo, parent_geo, fp_super_tile_geo, & - is_fine_pe, global_nest_domain, position, p_grid, n_grid, wt_h, istart_coarse, jstart_coarse) + is_fine_pe, global_nest_domain, position, p_grid, n_grid, wt_h, istart_coarse, jstart_coarse, Atm(child_grid_num)%neststruct%ind_h) call mn_meta_recalc( delta_i_c, delta_j_c, x_refine, y_refine, tile_geo_u, parent_geo, fp_super_tile_geo, & - is_fine_pe, global_nest_domain, position_u, p_grid_u, n_grid_u, wt_u, istart_coarse, jstart_coarse) + is_fine_pe, global_nest_domain, position_u, p_grid_u, n_grid_u, wt_u, istart_coarse, jstart_coarse, Atm(child_grid_num)%neststruct%ind_u) call mn_meta_recalc( delta_i_c, delta_j_c, x_refine, y_refine, tile_geo_v, parent_geo, fp_super_tile_geo, & - is_fine_pe, global_nest_domain, position_v, p_grid_v, n_grid_v, wt_v, istart_coarse, jstart_coarse) + is_fine_pe, global_nest_domain, position_v, p_grid_v, n_grid_v, wt_v, istart_coarse, jstart_coarse, Atm(child_grid_num)%neststruct%ind_v) if (use_timers) call mpp_clock_end (id_movnest5_3) endif @@ -936,10 +950,10 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, !! Step 5.3 -- Adjust the indices by the values of delta_i_c, delta_j_c !!============================================================================ - call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_h) - call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_u) - call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_v) - call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_b) + !call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_h) + !call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_u) + !call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_v) + !call mn_shift_index(delta_i_c, delta_j_c, Atm(child_grid_num)%neststruct%ind_b) if (debug_sync) call mpp_sync(full_pelist) ! Used to make debugging easier. Can be removed. @@ -996,6 +1010,8 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, case (2) ! Static nest smoothing algorithm - interpolation of coarse terrain in halo zone and 5 point blending zone of coarse and fine data call set_blended_terrain(Atm(n), mn_static%parent_orog_grid, mn_static%orog_grid, x_refine, Atm(n)%bd%ng, 10, a_step) + case (4) ! Use coarse terrain; no-op here. + ; case (5) ! 5 pt smoother. blend zone of 5 to match static nest call set_smooth_nest_terrain(Atm(n), mn_static%orog_grid, x_refine, 5, Atm(n)%bd%ng, 5) diff --git a/moving_nest/fv_moving_nest_types.F90 b/moving_nest/fv_moving_nest_types.F90 index 0a2a088fc..45bd12504 100644 --- a/moving_nest/fv_moving_nest_types.F90 +++ b/moving_nest/fv_moving_nest_types.F90 @@ -46,7 +46,7 @@ module fv_moving_nest_types_mod ! Moving Nest Namelist Variables logical :: is_moving_nest = .false. character(len=120) :: surface_dir = "INPUT/moving_nest" - integer :: terrain_smoother = 1 + integer :: terrain_smoother = 4 integer :: vortex_tracker = 0 integer :: ntrack = 1 integer :: corral_x = 5 @@ -221,7 +221,7 @@ module fv_moving_nest_types_mod ! Moving Nest Namelist Variables logical, dimension(MAX_NNEST) :: is_moving_nest = .False. character(len=120) :: surface_dir = "INPUT/moving_nest" - integer, dimension(MAX_NNEST) :: terrain_smoother = 1 ! 0 -- all high-resolution data, 1 - static nest smoothing algorithm with blending zone of 5 points, 2 - blending zone of 10 points, 5 - 5 point smoother, 9 - 9 point smoother + integer, dimension(MAX_NNEST) :: terrain_smoother = 4 ! 0 -- all high-resolution data, 1 - static nest smoothing algorithm with blending zone of 5 points, 2 - blending zone of 10 points, 5 - 5 point smoother, 9 - 9 point smoother integer, dimension(MAX_NNEST) :: vortex_tracker = 0 ! 0 - not a moving nest, tracker not needed ! 1 - prescribed nest moving ! 2 - following child domain center diff --git a/moving_nest/fv_moving_nest_utils.F90 b/moving_nest/fv_moving_nest_utils.F90 index 609a327ce..11485aa88 100644 --- a/moving_nest/fv_moving_nest_utils.F90 +++ b/moving_nest/fv_moving_nest_utils.F90 @@ -930,7 +930,7 @@ subroutine load_nest_latlons_from_nc(nc_filename, nxp, nyp, refine, pelist, & character(*), intent(in) :: nc_filename integer, intent(in) :: nxp, nyp, refine integer, allocatable, intent(in) :: pelist(:) - type(grid_geometry), intent(out) :: fp_tile_geo + type(grid_geometry), intent(inout) :: fp_tile_geo integer, intent(out) :: fp_istart_fine, fp_iend_fine, fp_jstart_fine, fp_jend_fine !======================================================================================== diff --git a/upp b/upp index e22724738..2b2c84a60 160000 --- a/upp +++ b/upp @@ -1 +1 @@ -Subproject commit e22724738fd104327fee7c3c7ffc805ccabd619f +Subproject commit 2b2c84a609f575fc293a4f90238391681bc383f0