From c5ddec2a68d95228cab8895b5dac1432db18d006 Mon Sep 17 00:00:00 2001 From: William Ramstrom Date: Fri, 15 Apr 2022 21:13:21 +0000 Subject: [PATCH 01/29] Performance optimization of moving nest. --- moving_nest/fv_moving_nest.F90 | 148 +++++++++++++++++-------- moving_nest/fv_moving_nest_main.F90 | 12 +- moving_nest/fv_moving_nest_physics.F90 | 10 +- moving_nest/fv_moving_nest_utils.F90 | 15 ++- 4 files changed, 129 insertions(+), 56 deletions(-) diff --git a/moving_nest/fv_moving_nest.F90 b/moving_nest/fv_moving_nest.F90 index 8fc645b6d..2d094132c 100644 --- a/moving_nest/fv_moving_nest.F90 +++ b/moving_nest/fv_moving_nest.F90 @@ -591,7 +591,69 @@ subroutine mn_var_fill_intern_nest_halos_r8_4d(data_var, domain_fine, is_fine_pe end subroutine mn_var_fill_intern_nest_halos_r8_4d + !>@brief Find the parent point that corresponds to the is,js point of the nest, and returns that nest point also + subroutine calc_nest_alignment(Atm, n, nest_x, nest_y, parent_x, parent_y) + type(fv_atmos_type), allocatable, intent(in) :: Atm(:) !< Atm data array + integer, intent(in) :: n !< Grid numbers + integer, intent(out) :: nest_x, nest_y, parent_x, parent_y + + integer :: refine + integer :: child_grid_num + integer :: ioffset, joffset + + child_grid_num = n + + refine = Atm(child_grid_num)%neststruct%refinement + + ! parent_x and parent_y are on the supergrid, so an increment of ioffset is an increment of 2*refine + + nest_x = Atm(child_grid_num)%bd%isd + nest_y = Atm(child_grid_num)%bd%jsd + + ioffset = Atm(n)%neststruct%ioffset + joffset = Atm(n)%neststruct%joffset + + ! Increment of 3 is for halo. Factor of 2 is for supergrid. + parent_x = (nest_x - 3)*2 + ioffset*refine*2 + parent_y = (nest_y - 3)*2 + joffset*refine*2 + + end subroutine calc_nest_alignment + + + + subroutine check_nest_alignment(nest_geo, parent_geo, nest_x, nest_y, parent_x, parent_y, found) + type(grid_geometry), intent(in) :: nest_geo !< Tile geometry + type(grid_geometry), intent(in) :: parent_geo !< Parent grid at high-resolution geometry + integer, intent(in) :: nest_x, nest_y, parent_x, parent_y + logical, intent(out) :: found + + real(kind=R_GRID) :: pi = 4 * atan(1.0d0) + real :: rad2deg + integer :: this_pe + + this_pe = mpp_pe() + rad2deg = 180.0 / pi + + found = .False. + if (abs(parent_geo%lats(parent_x, parent_y) - nest_geo%lats(nest_x, nest_y)) .lt. 0.0001) then + if (abs(parent_geo%lons(parent_x, parent_y) - nest_geo%lons(nest_x, nest_y)) .lt. 0.0001) then + found = .True. + endif + if (abs(abs(parent_geo%lons(parent_x, parent_y) - nest_geo%lons(nest_x, nest_y)) - 2*pi) .lt. 0.0001) then + found = .True. + endif + endif + + !print '("[INFO] WDR C-ALIGN check_nest_alignment npe=",I0," found=",L1," parent(",I0,",",I0,") nest(",I0,",",I0,")",4F16.9)', & + ! this_pe, found, parent_x, parent_y, nest_x, nest_y, & + ! parent_geo%lats(parent_x, parent_y)*rad2deg, parent_geo%lons(parent_x, parent_y)*rad2deg, & + ! nest_geo%lats(nest_x, nest_y)*rad2deg, nest_geo%lons(nest_x, nest_y)*rad2deg + + + + end subroutine check_nest_alignment + !!============================================================================ !! Step 5.1 -- Load the latlon data from NetCDF !! update parent_geo, tile_geo*, p_grid*, n_grid* @@ -689,12 +751,6 @@ subroutine mn_latlon_load_parent(surface_dir, Atm, n, parent_tile, delta_i_c, de enddo if (debug_log) call show_tile_geo(tile_geo, this_pe, "tile_geo") - call find_nest_alignment(tile_geo, fp_super_tile_geo, nest_x, nest_y, parent_x, parent_y) - - if (parent_x .eq. -999) then - print '("[ERROR] WDR mn_latlon_load_parent on npe=",I0," parent and nest grids are not aligned!")', this_pe - call mpp_error(FATAL, "mn_latlon_load_parent parent and nest grids are not aligned.") - endif ! Allocate tile_geo_u just for this PE, copied from Atm(n)%gridstruct%grid ! grid is 1 larger than agrid @@ -704,26 +760,18 @@ subroutine mn_latlon_load_parent(surface_dir, Atm, n, parent_tile, delta_i_c, de tile_geo_u%nxp = tile_geo_u%nx + 1 tile_geo_u%nyp = tile_geo_u%ny + 1 - allocate(tile_geo_u%lons(lbound(Atm(n)%gridstruct%agrid, 1):ubound(Atm(n)%gridstruct%agrid, 1), lbound(Atm(n)%gridstruct%grid, 2):ubound(Atm(n)%gridstruct%grid, 2))) - allocate(tile_geo_u%lats(lbound(Atm(n)%gridstruct%agrid, 1):ubound(Atm(n)%gridstruct%agrid, 1), lbound(Atm(n)%gridstruct%grid, 2):ubound(Atm(n)%gridstruct%grid, 2))) + + if (.not. allocated(tile_geo_u%lons)) then + !print '("[INFO] WDR mn_latlon_load_parent ALLOCATE tile_geo_u%lons npe=",I0)', this_pe + allocate(tile_geo_u%lons(lbound(Atm(n)%gridstruct%agrid, 1):ubound(Atm(n)%gridstruct%agrid, 1), lbound(Atm(n)%gridstruct%grid, 2):ubound(Atm(n)%gridstruct%grid, 2))) + allocate(tile_geo_u%lats(lbound(Atm(n)%gridstruct%agrid, 1):ubound(Atm(n)%gridstruct%agrid, 1), lbound(Atm(n)%gridstruct%grid, 2):ubound(Atm(n)%gridstruct%grid, 2))) + !else + ! print '("[INFO] WDR mn_latlon_load_parent ALREADY ALLOCATED tile_geo_u%lons npe=",I0)', this_pe + endif tile_geo_u%lons = -999.9 tile_geo_u%lats = -999.9 - do x = lbound(tile_geo_u%lats, 1), ubound(tile_geo_u%lats, 1) - do y = lbound(tile_geo_u%lats, 2), ubound(tile_geo_u%lats, 2) - fp_i = (x - nest_x) * 2 + parent_x - 1 - fp_j = (y - nest_y) * 2 + parent_y - - !print '("[INFO] WDR mn_latlon_load_parent on npe=",I0," fp_i=",I0," fp_j=",I0,4I6)', this_pe, fp_i, fp_j, nest_x, nest_y, parent_x, parent_y - - tile_geo_u%lons(x,y) = fp_super_tile_geo%lons(fp_i, fp_j) - tile_geo_u%lats(x,y) = fp_super_tile_geo%lats(fp_i, fp_j) - enddo - enddo - - if (debug_log) call show_tile_geo(tile_geo_u, this_pe, "tile_geo_u") - ! Allocate tile_geo_v just for this PE, copied from Atm(n)%gridstruct%grid ! grid is 1 larger than agrid ! u(npx, npy+1) @@ -738,18 +786,6 @@ subroutine mn_latlon_load_parent(surface_dir, Atm, n, parent_tile, delta_i_c, de tile_geo_v%lons = -999.9 tile_geo_v%lats = -999.9 - do x = lbound(tile_geo_v%lats, 1), ubound(tile_geo_v%lats, 1) - do y = lbound(tile_geo_v%lats, 2), ubound(tile_geo_v%lats, 2) - fp_i = (x - nest_x) * 2 + parent_x - fp_j = (y - nest_y) * 2 + parent_y - 1 - - tile_geo_v%lons(x,y) = fp_super_tile_geo%lons(fp_i, fp_j) - tile_geo_v%lats(x,y) = fp_super_tile_geo%lats(fp_i, fp_j) - enddo - enddo - - if (debug_log) call show_tile_geo(tile_geo_v, this_pe, "tile_geo_v") - !=========================================================== ! End tile_geo per PE. !=========================================================== @@ -775,7 +811,12 @@ subroutine mn_latlon_load_parent(surface_dir, Atm, n, parent_tile, delta_i_c, de if (debug_log) print '("[INFO] WDR MV_NST2 bounds1 (tile_geo%lats)=",I0,"-",I0)', lbound(tile_geo%lats,1), ubound(tile_geo%lats,1) if (debug_log) print '("[INFO] WDR MV_NST2 bounds2 (tile_geo%lats)=",I0,"-",I0)', lbound(tile_geo%lats,2), ubound(tile_geo%lats,2) - call move_nest_geo(tile_geo, tile_geo_u, tile_geo_v, fp_super_tile_geo, delta_i_c, delta_j_c, x_refine, y_refine) + call move_nest_geo(Atm, n, tile_geo, tile_geo_u, tile_geo_v, fp_super_tile_geo, delta_i_c, delta_j_c, x_refine, y_refine) + + if (parent_x .eq. -999) then + print '("[ERROR] WDR mn_latlon_load_parent on npe=",I0," parent and nest grids are not aligned!")', this_pe + call mpp_error(FATAL, "mn_latlon_load_parent parent and nest grids are not aligned.") + endif call assign_n_p_grids(parent_geo, tile_geo, p_grid, n_grid, position) call assign_n_p_grids(parent_geo, tile_geo_u, p_grid_u, n_grid_u, position_u) @@ -1682,7 +1723,7 @@ subroutine mn_meta_reset_gridstruct(Atm, n, child_grid_num, nest_domain, fp_supe logical, save :: first_time = .true. integer, save :: id_reset1, id_reset2, id_reset3, id_reset4, id_reset5, id_reset6, id_reset7 - logical :: use_timers = .false. ! Set this to true to generate performance profiling information in out.* file + logical :: use_timers = .False. ! Set this to true to generate performance profiling information in out.* file if (first_time .and. use_timers) then id_reset1 = mpp_clock_id ('MN 7 Reset 1', flags = clock_flag_default, grain=CLOCK_ROUTINE ) @@ -2478,33 +2519,43 @@ end function almost_equal !>@brief The subroutine 'move_nest_geo' shifts tile_geo values using the data from fp_super_tile_geo - subroutine move_nest_geo(tile_geo, tile_geo_u, tile_geo_v, fp_super_tile_geo, delta_i_c, delta_j_c, x_refine, y_refine) + subroutine move_nest_geo(Atm, n, tile_geo, tile_geo_u, tile_geo_v, fp_super_tile_geo, delta_i_c, delta_j_c, x_refine, y_refine) implicit none + type(fv_atmos_type), allocatable, intent(in) :: Atm(:) !< Atm data array + integer, intent(in) :: n !< Grid numbers type(grid_geometry), intent(inout) :: tile_geo !< A-grid tile geometry type(grid_geometry), intent(inout) :: tile_geo_u !< u-wind tile geometry type(grid_geometry), intent(inout) :: tile_geo_v !< v-wind tile geometry type(grid_geometry), intent(in) :: fp_super_tile_geo !< Parent high-resolution supergrid tile geometry integer, intent(in) :: delta_i_c, delta_j_c, x_refine, y_refine !< delta i,j for nest move. Nest refinement. - integer :: nest_x, nest_y, parent_x, parent_y - - type(bbox) :: tile_bbox, fp_tile_bbox, tile_bbox_u, tile_bbox_v - integer :: i, j, fp_i, fp_j + integer :: nest_x, nest_y, parent_x, parent_y + type(bbox) :: tile_bbox, fp_tile_bbox, tile_bbox_u, tile_bbox_v + integer :: i, j, fp_i, fp_j + integer :: this_pe + logical :: found ! tile_geo is cell-centered, at nest refinement ! fp_super_tile_geo is a supergrid, at nest refinement - call find_nest_alignment(tile_geo, fp_super_tile_geo, nest_x, nest_y, parent_x, parent_y) + this_pe = mpp_pe() call fill_bbox(tile_bbox, tile_geo%lats) call fill_bbox(tile_bbox_u, tile_geo_u%lats) call fill_bbox(tile_bbox_v, tile_geo_v%lats) call fill_bbox(fp_tile_bbox, fp_super_tile_geo%lats) - ! Calculate new parent alignment -- supergrid at the refine ratio - ! delta_{i,j}_c are at the coarse center grid resolution - parent_x = parent_x + delta_i_c * 2 * x_refine - parent_y = parent_y + delta_j_c * 2 * y_refine + !! Calculate new parent alignment -- supergrid at the refine ratio + !! delta_{i,j}_c are at the coarse center grid resolution + !parent_x = parent_x + delta_i_c * 2 * x_refine + !parent_y = parent_y + delta_j_c * 2 * y_refine + + !print '("[INFO] WDR ALIGN-D npe=",I0," ioffset=",I0," joffset=",I0," delta_i_c=",I0," delta_j_c=",I0," nest_x=",I0," nest_y=",I0," parent_x=",I0," parent_y=",I0)', this_pe, ioffset, joffset, delta_i_c, delta_j_c, nest_x, nest_y, parent_x, parent_y + + call calc_nest_alignment(Atm, n, nest_x, nest_y, parent_x, parent_y) + + !print '("[INFO] WDR ALIGN-E npe=",I0," ioffset=",I0," joffset=",I0," delta_i_c=",I0," delta_j_c=",I0," nest_x=",I0," nest_y=",I0," parent_x=",I0," parent_y=",I0)', & + ! this_pe, Atm(n)%neststruct%ioffset, Atm(n)%neststruct%joffset, delta_i_c, delta_j_c, nest_x, nest_y, parent_x, parent_y ! Brute force repopulation of full tile_geo grids. ! Optimization would be to use EOSHIFT and bring in just leading edge @@ -2566,7 +2617,10 @@ subroutine move_nest_geo(tile_geo, tile_geo_u, tile_geo_v, fp_super_tile_geo, de enddo ! Validate at the end - call find_nest_alignment(tile_geo, fp_super_tile_geo, nest_x, nest_y, parent_x, parent_y) + call check_nest_alignment(tile_geo, fp_super_tile_geo, nest_x, nest_y, parent_x, parent_y, found) + + !print '("[INFO] WDR ALIGN-C npe=",I0," delta_i_c=",I0," delta_j_c=",I0," nest_x=",I0," nest_y=",I0," parent_x=",I0," parent_y=",I0)', this_pe, delta_i_c, delta_j_c, nest_x, nest_y, parent_x, parent_y + end subroutine move_nest_geo diff --git a/moving_nest/fv_moving_nest_main.F90 b/moving_nest/fv_moving_nest_main.F90 index c8691db18..2cdd1ac83 100644 --- a/moving_nest/fv_moving_nest_main.F90 +++ b/moving_nest/fv_moving_nest_main.F90 @@ -125,7 +125,7 @@ module fv_moving_nest_main_mod use fv_moving_nest_physics_mod, only: mn_reset_phys_latlon, mn_surface_grids ! Grid reset routines - use fv_moving_nest_mod, only: grid_geometry, assign_n_p_grids, move_nest_geo + use fv_moving_nest_mod, only: grid_geometry use fv_moving_nest_utils_mod, only: fill_grid_from_supergrid, fill_weight_grid ! Physics moving logical variables @@ -494,7 +494,7 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, real(kind=R_GRID), allocatable :: p_grid(:,:,:), n_grid(:,:,:) real(kind=R_GRID), allocatable :: p_grid_u(:,:,:), n_grid_u(:,:,:) real(kind=R_GRID), allocatable :: p_grid_v(:,:,:), n_grid_v(:,:,:) - real, allocatable :: wt_h(:,:,:) + real, allocatable :: wt_h(:,:,:) ! TODO verify that these are deallocated real, allocatable :: wt_u(:,:,:) real, allocatable :: wt_v(:,:,:) !real :: ua(isd:ied,jsd:jed) @@ -1174,6 +1174,14 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, !call compare_terrain("phis", Atm(n)%phis, 1, Atm(n)%neststruct%ind_h, x_refine, y_refine, is_fine_pe, global_nest_domain) + !deallocate(tile_geo%lats, tile_geo%lons) + !deallocate(tile_geo_u%lats, tile_geo_u%lons) + !deallocate(tile_geo_v%lats, tile_geo_v%lons) + + !deallocate(p_grid, n_grid) + !deallocate(p_grid_u, n_grid_u) + !deallocate(p_grid_v, n_grid_v) + if (debug_log) call show_nest_grid(Atm(n), this_pe, 99) end subroutine fv_moving_nest_exec diff --git a/moving_nest/fv_moving_nest_physics.F90 b/moving_nest/fv_moving_nest_physics.F90 index 074bb3c16..c2ff732fc 100644 --- a/moving_nest/fv_moving_nest_physics.F90 +++ b/moving_nest/fv_moving_nest_physics.F90 @@ -94,7 +94,7 @@ module fv_moving_nest_physics_mod use fv_moving_nest_utils_mod, only: fill_nest_from_buffer_cell_center_masked use fv_moving_nest_utils_mod, only: fill_nest_halos_from_parent_masked - use fv_moving_nest_mod, only: mn_var_fill_intern_nest_halos, mn_var_dump_to_netcdf, mn_var_shift_data + use fv_moving_nest_mod, only: mn_var_fill_intern_nest_halos, mn_var_dump_to_netcdf, mn_var_shift_data, calc_nest_alignment use fv_moving_nest_types_mod, only: Moving_nest implicit none @@ -227,7 +227,7 @@ end subroutine mn_phys_reset_sfc_props !>@brief The subroutine 'mn_phys_reset_phys_latlon' sets the lat/lons from the high-resolution input file data !>@details This subroutine sets lat/lons of the moved nest, then recalculates all the derived quantities (dx,dy,etc.) subroutine mn_reset_phys_latlon(Atm, n, tile_geo, fp_super_tile_geo, Atm_block, IPD_control, IPD_data) - type(fv_atmos_type), intent(in) :: Atm(:) !< Array of atmospheric data + type(fv_atmos_type), allocatable, intent(in) :: Atm(:) !< Array of atmospheric data integer, intent(in) :: n !< Current grid number type(grid_geometry), intent(in) :: tile_geo !< Bounds of this grid type(grid_geometry), intent(in) :: fp_super_tile_geo !< Bounds of high-resolution parent grid @@ -253,10 +253,8 @@ subroutine mn_reset_phys_latlon(Atm, n, tile_geo, fp_super_tile_geo, Atm_block, allocate(lons(isc:iec, jsc:jec)) allocate(area(isc:iec, jsc:jec)) - ! This is going to be slow -- replace with better way to calculate index offsets, or pass them from earlier calculations - ! TODO optimization here - call find_nest_alignment(tile_geo, fp_super_tile_geo, nest_x, nest_y, parent_x, parent_y) - !print '("WDR mn_reset_phys_latlon AB npe=",I0)', this_pe + call calc_nest_alignment(Atm, n, nest_x, nest_y, parent_x, parent_y) + !print '("[INFO] WDR ALIGN-PHYS-NEW npe=",I0," nest_x=",I0," nest_y=",I0," parent_x=",I0," parent_y=",I0)', this_pe, nest_x, nest_y, parent_x, parent_y do x = isc, iec do y = jsc, jec diff --git a/moving_nest/fv_moving_nest_utils.F90 b/moving_nest/fv_moving_nest_utils.F90 index d751214cc..aa64875a8 100644 --- a/moving_nest/fv_moving_nest_utils.F90 +++ b/moving_nest/fv_moving_nest_utils.F90 @@ -30,7 +30,7 @@ module fv_moving_nest_utils_mod #ifdef MOVING_NEST - + use fms_mod, only : mpp_clock_id, mpp_clock_begin, mpp_clock_end, CLOCK_ROUTINE, clock_flag_default use mpp_mod, only: FATAL, WARNING, MPP_DEBUG, NOTE, MPP_CLOCK_SYNC,MPP_CLOCK_DETAILED use mpp_mod, only: mpp_pe, mpp_npes, mpp_root_pe, mpp_error, mpp_set_warn_level use mpp_mod, only: mpp_declare_pelist, mpp_set_current_pelist, mpp_sync, mpp_sync_self @@ -1582,8 +1582,18 @@ subroutine find_nest_alignment(nest_geo, parent_geo, nest_x, nest_y, parent_x, p real(kind=R_GRID) :: pi = 4 * atan(1.0d0) real :: rad2deg integer :: this_pe + + logical, save :: first_time = .true. + integer, save :: id_nest_align + this_pe = mpp_pe() + + if (first_time) then + id_nest_align = mpp_clock_id ('MN Nest Align', flags = clock_flag_default, grain=CLOCK_ROUTINE ) + first_time = .false. + endif + call mpp_clock_begin (id_nest_align) rad2deg = 180.0 / pi @@ -1655,6 +1665,9 @@ subroutine find_nest_alignment(nest_geo, parent_geo, nest_x, nest_y, parent_x, p enddo endif + call mpp_clock_end (id_nest_align) + + end subroutine find_nest_alignment From 2b413b549d344522cd6eb4db0c53d8d68081f1cd Mon Sep 17 00:00:00 2001 From: William Ramstrom Date: Mon, 2 May 2022 02:35:18 +0000 Subject: [PATCH 02/29] Moving nest performance optimization stage 2. --- moving_nest/fv_moving_nest.F90 | 177 ++++++++++++++++++++----- moving_nest/fv_moving_nest_main.F90 | 33 ++++- moving_nest/fv_moving_nest_physics.F90 | 1 - tools/fv_grid_tools.F90 | 13 +- 4 files changed, 177 insertions(+), 47 deletions(-) diff --git a/moving_nest/fv_moving_nest.F90 b/moving_nest/fv_moving_nest.F90 index 2d094132c..1e1f2a8be 100644 --- a/moving_nest/fv_moving_nest.F90 +++ b/moving_nest/fv_moving_nest.F90 @@ -56,7 +56,7 @@ module fv_moving_nest_mod #ifdef MOVING_NEST 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 + 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_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 @@ -668,9 +668,12 @@ subroutine mn_latlon_load_parent(surface_dir, Atm, n, parent_tile, delta_i_c, de integer, intent(in) :: delta_i_c, delta_j_c !< Nest motion in delta i,j type(grid_geometry), intent(inout) :: parent_geo, tile_geo, tile_geo_u, tile_geo_v !< Tile geometries type(grid_geometry), intent(in) :: fp_super_tile_geo !< Parent grid at high-resolution geometry - real(kind=R_GRID), allocatable, intent(out) :: p_grid(:,:,:), n_grid(:,:,:) !< A-stagger lat/lon grids - real(kind=R_GRID), allocatable, intent(out) :: p_grid_u(:,:,:), n_grid_u(:,:,:) !< u-wind staggered lat/lon grids - real(kind=R_GRID), allocatable, intent(out) :: p_grid_v(:,:,:), n_grid_v(:,:,:) !< v-wind staggered lat/lon grids + real(kind=R_GRID), allocatable, intent(inout):: p_grid(:,:,:) !< A-stagger lat/lon grids + real(kind=R_GRID), allocatable, intent(inout):: p_grid_u(:,:,:) !< u-wind staggered lat/lon grids + real(kind=R_GRID), allocatable, intent(inout):: p_grid_v(:,:,:) !< v-wind staggered lat/lon grids + real(kind=R_GRID), allocatable, intent(out) :: n_grid(:,:,:) !< A-stagger lat/lon grids + real(kind=R_GRID), allocatable, intent(out) :: n_grid_u(:,:,:) !< u-wind staggered lat/lon grids + real(kind=R_GRID), allocatable, intent(out) :: n_grid_v(:,:,:) !< v-wind staggered lat/lon grids character(len=256) :: grid_filename logical, save :: first_nest_move = .true. @@ -681,8 +684,22 @@ subroutine mn_latlon_load_parent(surface_dir, Atm, n, parent_tile, delta_i_c, de integer :: nest_x, nest_y, parent_x, parent_y integer :: this_pe + logical, save :: first_time = .True. + integer, save :: id_load1, id_load2, id_load3, id_load4, id_load5 + logical :: use_timers = .True. + this_pe = mpp_pe() + if (first_time) then + id_load1 = mpp_clock_id ('MN LatLon Part 1 File', flags = clock_flag_default, grain=CLOCK_SUBCOMPONENT ) + id_load2 = mpp_clock_id ('MN LatLon Part 2 File', flags = clock_flag_default, grain=CLOCK_SUBCOMPONENT ) + id_load3 = mpp_clock_id ('MN LatLon Part 3 File', flags = clock_flag_default, grain=CLOCK_SUBCOMPONENT ) + id_load4 = mpp_clock_id ('MN LatLon Part 4 File', flags = clock_flag_default, grain=CLOCK_SUBCOMPONENT ) + id_load5 = mpp_clock_id ('MN LatLon Part 5 File', flags = clock_flag_default, grain=CLOCK_SUBCOMPONENT ) + + first_time = .False. + endif + position = CENTER position_u = NORTH position_v = EAST @@ -696,14 +713,28 @@ subroutine mn_latlon_load_parent(surface_dir, Atm, n, parent_tile, delta_i_c, de if (first_nest_move) then if (debug_log) print '("[INFO] WDR mn_latlon_load_parent READING static coarse file on npe=",I0)', this_pe + if (use_timers) call mpp_clock_begin (id_load1) 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, & 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)) + + ! These are big (parent grid size), and do not change during the model integration. + call assign_p_grids(parent_geo, p_grid, position) + call assign_p_grids(parent_geo, p_grid_u, position_u) + call assign_p_grids(parent_geo, p_grid_v, position_v) + first_nest_move = .false. + if (use_timers) call mpp_clock_end (id_load1) endif + if (use_timers) call mpp_clock_begin (id_load2) + parent_geo%nxp = Atm(1)%npx parent_geo%nyp = Atm(1)%npy @@ -750,6 +781,10 @@ subroutine mn_latlon_load_parent(surface_dir, Atm, n, parent_tile, delta_i_c, de enddo enddo + if (use_timers) call mpp_clock_end (id_load2) + if (use_timers) call mpp_clock_begin (id_load3) + + if (debug_log) call show_tile_geo(tile_geo, this_pe, "tile_geo") ! Allocate tile_geo_u just for this PE, copied from Atm(n)%gridstruct%grid @@ -790,15 +825,12 @@ subroutine mn_latlon_load_parent(surface_dir, Atm, n, parent_tile, delta_i_c, de ! End tile_geo per PE. !=========================================================== - allocate(p_grid(1:parent_geo%nxp, 1:parent_geo%nyp,2)) allocate(n_grid(Atm(child_grid_num)%bd%isd:Atm(child_grid_num)%bd%ied, Atm(child_grid_num)%bd%jsd:Atm(child_grid_num)%bd%jed, 2)) n_grid = real_snan - allocate(p_grid_u(1:parent_geo%nxp, 1:parent_geo%nyp+1,2)) allocate(n_grid_u(Atm(child_grid_num)%bd%isd:Atm(child_grid_num)%bd%ied, Atm(child_grid_num)%bd%jsd:Atm(child_grid_num)%bd%jed+1, 2)) n_grid_u = real_snan - allocate(p_grid_v(1:parent_geo%nxp+1, 1:parent_geo%nyp,2)) allocate(n_grid_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, 2)) n_grid_v = real_snan @@ -811,16 +843,27 @@ subroutine mn_latlon_load_parent(surface_dir, Atm, n, parent_tile, delta_i_c, de if (debug_log) print '("[INFO] WDR MV_NST2 bounds1 (tile_geo%lats)=",I0,"-",I0)', lbound(tile_geo%lats,1), ubound(tile_geo%lats,1) if (debug_log) print '("[INFO] WDR MV_NST2 bounds2 (tile_geo%lats)=",I0,"-",I0)', lbound(tile_geo%lats,2), ubound(tile_geo%lats,2) + + if (use_timers) call mpp_clock_end (id_load3) + if (use_timers) call mpp_clock_begin (id_load4) + call move_nest_geo(Atm, n, tile_geo, tile_geo_u, tile_geo_v, fp_super_tile_geo, delta_i_c, delta_j_c, x_refine, y_refine) + if (use_timers) call mpp_clock_end (id_load4) + if (use_timers) call mpp_clock_begin (id_load5) + + if (parent_x .eq. -999) then print '("[ERROR] WDR mn_latlon_load_parent on npe=",I0," parent and nest grids are not aligned!")', this_pe call mpp_error(FATAL, "mn_latlon_load_parent parent and nest grids are not aligned.") endif - call assign_n_p_grids(parent_geo, tile_geo, p_grid, n_grid, position) - call assign_n_p_grids(parent_geo, tile_geo_u, p_grid_u, n_grid_u, position_u) - call assign_n_p_grids(parent_geo, tile_geo_v, p_grid_v, n_grid_v, position_v) + ! These grids are small (nest size), and change each time nest moves. + call assign_n_grids(tile_geo, n_grid, position) + call assign_n_grids(tile_geo_u, n_grid_u, position_u) + call assign_n_grids(tile_geo_v, n_grid_v, position_v) + + if (use_timers) call mpp_clock_end (id_load5) end subroutine mn_latlon_load_parent @@ -1187,14 +1230,10 @@ subroutine mn_var_shift_data_r4_2d(data_var, interp_type, wt, ind, delta_i_c, de integer, intent(in) :: position !< Grid offset real*4, dimension(:,:), allocatable :: nbuffer, sbuffer, ebuffer, wbuffer - logical :: parent_proc, child_proc type(bbox) :: north_fine, north_coarse ! step 4 type(bbox) :: south_fine, south_coarse type(bbox) :: east_fine, east_coarse type(bbox) :: west_fine, west_coarse - integer :: my_stat - character(256) :: my_errmsg - integer :: is, ie, js, je integer :: this_pe integer :: nest_level = 1 ! WDR TODO allow to vary @@ -1302,14 +1341,10 @@ subroutine mn_var_shift_data_r8_2d(data_var, interp_type, wt, ind, delta_i_c, de integer, intent(in) :: position !< Grid offset real*8, dimension(:,:), allocatable :: nbuffer, sbuffer, ebuffer, wbuffer - logical :: parent_proc, child_proc type(bbox) :: north_fine, north_coarse ! step 4 type(bbox) :: south_fine, south_coarse type(bbox) :: east_fine, east_coarse type(bbox) :: west_fine, west_coarse - integer :: my_stat - character(256) :: my_errmsg - integer :: is, ie, js, je integer :: this_pe integer :: nest_level = 1 ! WDR TODO allow to vary @@ -1382,14 +1417,10 @@ subroutine mn_var_shift_data_r4_3d(data_var, interp_type, wt, ind, delta_i_c, de integer, intent(in) :: position, nz !< Grid offset, number of vertical levels real*4, dimension(:,:,:), allocatable :: nbuffer, sbuffer, ebuffer, wbuffer - logical :: parent_proc, child_proc type(bbox) :: north_fine, north_coarse ! step 4 type(bbox) :: south_fine, south_coarse type(bbox) :: east_fine, east_coarse type(bbox) :: west_fine, west_coarse - integer :: my_stat - character(256) :: my_errmsg - integer :: is, ie, js, je integer :: this_pe integer :: nest_level = 1 ! WDR TODO allow to vary @@ -1462,14 +1493,10 @@ subroutine mn_var_shift_data_r8_3d(data_var, interp_type, wt, ind, delta_i_c, de integer, intent(in) :: position, nz !< Grid offset, number vertical levels real*8, dimension(:,:,:), allocatable :: nbuffer, sbuffer, ebuffer, wbuffer - logical :: parent_proc, child_proc type(bbox) :: north_fine, north_coarse ! step 4 type(bbox) :: south_fine, south_coarse type(bbox) :: east_fine, east_coarse type(bbox) :: west_fine, west_coarse - integer :: my_stat - character(256) :: my_errmsg - integer :: is, ie, js, je integer :: this_pe integer :: nest_level = 1 ! WDR TODO allow to vary @@ -1539,16 +1566,12 @@ subroutine mn_var_shift_data_r4_4d(data_var, interp_type, wt, ind, delta_i_c, de integer, intent(in) :: position, nz !< Grid offset, number of vertical levels real*4, dimension(:,:,:,:), allocatable :: nbuffer, sbuffer, ebuffer, wbuffer - logical :: parent_proc, child_proc type(bbox) :: north_fine, north_coarse ! step 4 type(bbox) :: south_fine, south_coarse type(bbox) :: east_fine, east_coarse type(bbox) :: west_fine, west_coarse - integer :: my_stat - character(256) :: my_errmsg integer :: n4d integer :: this_pe - integer :: is, ie, js, je integer :: nest_level = 1 ! WDR TODO allow to vary this_pe = mpp_pe() @@ -1619,16 +1642,12 @@ subroutine mn_var_shift_data_r8_4d(data_var, interp_type, wt, ind, delta_i_c, de integer, intent(in) :: position, nz !< Grid offset, number of vertical levels real*8, dimension(:,:,:,:), allocatable :: nbuffer, sbuffer, ebuffer, wbuffer - logical :: parent_proc, child_proc type(bbox) :: north_fine, north_coarse ! step 4 type(bbox) :: south_fine, south_coarse type(bbox) :: east_fine, east_coarse type(bbox) :: west_fine, west_coarse - integer :: my_stat - character(256) :: my_errmsg integer :: n4d integer :: this_pe - integer :: is, ie, js, je integer :: nest_level = 1 ! WDR TODO allow to vary this_pe = mpp_pe() @@ -1723,7 +1742,7 @@ subroutine mn_meta_reset_gridstruct(Atm, n, child_grid_num, nest_domain, fp_supe logical, save :: first_time = .true. integer, save :: id_reset1, id_reset2, id_reset3, id_reset4, id_reset5, id_reset6, id_reset7 - logical :: use_timers = .False. ! Set this to true to generate performance profiling information in out.* file + logical :: use_timers = .True. ! Set this to true to generate performance profiling information in out.* file if (first_time .and. use_timers) then id_reset1 = mpp_clock_id ('MN 7 Reset 1', flags = clock_flag_default, grain=CLOCK_ROUTINE ) @@ -2693,6 +2712,94 @@ subroutine assign_n_p_grids(parent_geo, tile_geo, p_grid, n_grid, position) end subroutine assign_n_p_grids + !>@brief The subroutine 'assign_p_grids' sets values for parent grid arrays from the grid_geometry structures. This is static through the model run. + subroutine assign_p_grids(parent_geo, p_grid, position) + type(grid_geometry), intent(in) :: parent_geo !< Parent geometry + real(kind=R_GRID), allocatable, intent(inout) :: p_grid(:,:,:) !< Parent grid + integer, intent(in) :: position !< Grid offset + + integer :: i,j + + if (position == CENTER) then + 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) + p_grid(i, j, 2) = parent_geo%lats(2*i, 2*j) + enddo + enddo + + ! 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) + p_grid(i, j, 2) = parent_geo%lats(2*i, 2*j-1) + enddo + enddo + + ! 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) + enddo + enddo + endif + + end subroutine assign_p_grids + + + + !>@brief The subroutine 'assign_n_grids' sets values for nest grid arrays from the grid_geometry structures. + subroutine assign_n_grids(tile_geo, n_grid, position) + type(grid_geometry), intent(in) :: tile_geo !< Parent geometry, nest geometry + real(kind=R_GRID), allocatable, intent(inout) :: n_grid(:,:,:) !< Nest grid + integer, intent(in) :: position !< Grid offset + + integer :: i,j + + if (position == CENTER) then + do j = lbound(tile_geo%lats,2), ubound(tile_geo%lats,2) + do i = lbound(tile_geo%lats,1), ubound(tile_geo%lats,1) + ! centered grid version + n_grid(i, j, 1) = tile_geo%lons(i, j) + n_grid(i, j, 2) = tile_geo%lats(i, j) + !if (debug_log) print '("[INFO] WDR populate ngrid npe=",I0, I4,I4, F12.4, F12.4)', this_pe, i, j, n_grid(i,j,1), n_grid(i,j,2) + enddo + enddo + + ! u(npx, npy+1) + elseif (position == NORTH) then ! u wind on D-stagger + do j = lbound(tile_geo%lats,2), ubound(tile_geo%lats,2) + do i = lbound(tile_geo%lats,1), ubound(tile_geo%lats,1) + ! centered grid version + n_grid(i, j, 1) = tile_geo%lons(i, j) + n_grid(i, j, 2) = tile_geo%lats(i, j) + !if (debug_log) print '("[INFO] WDR populate ngrid_u npe=",I0, I4,I4, F12.4, F12.4)', this_pe, i, j, n_grid(i,j,1), n_grid(i,j,2) + enddo + enddo + + ! v(npx+1, npy) + elseif (position == EAST) then ! v wind on D-stagger + do j = lbound(tile_geo%lats,2), ubound(tile_geo%lats,2) + do i = lbound(tile_geo%lats,1), ubound(tile_geo%lats,1) + ! centered grid version + n_grid(i, j, 1) = tile_geo%lons(i, j) + n_grid(i, j, 2) = tile_geo%lats(i, j) + !if (debug_log) print '("[INFO] WDR populate ngrid_v npe=",I0, I4,I4, F12.4, F12.4)', this_pe, i, j, n_grid(i,j,1), n_grid(i,j,2) + enddo + enddo + + endif + + end subroutine assign_n_grids + + + !>@brief The subroutine 'calc_nest_halo_weights' calculates the interpolation weights !>@details Computationally demanding; target for optimization after nest moves diff --git a/moving_nest/fv_moving_nest_main.F90 b/moving_nest/fv_moving_nest_main.F90 index 2cdd1ac83..1a78dc0ff 100644 --- a/moving_nest/fv_moving_nest_main.F90 +++ b/moving_nest/fv_moving_nest_main.F90 @@ -161,9 +161,10 @@ module fv_moving_nest_main_mod ! --- Clock ids for moving_nest performance metering integer :: id_movnest1, id_movnest1_9, id_movnest2, id_movnest3, id_movnest4, id_movnest5 + integer :: id_movnest5_1, id_movnest5_2, id_movnest5_3, id_movnest5_4 integer :: id_movnest6, id_movnest7_0, id_movnest7_1, id_movnest7_2, id_movnest7_3, id_movnest8, id_movnest9 integer :: id_movnestTot - logical :: use_timers = .False. ! Set this to true for detailed performance profiling. False only profiles total moving nest time. + logical :: use_timers = .True. ! Set this to true for detailed performance profiling. False only profiles total moving nest time. integer, save :: output_step = 0 contains @@ -260,6 +261,11 @@ subroutine fv_moving_nest_init_clocks() id_movnest3 = mpp_clock_id ('MN Part 3 Meta Move Nest', flags = clock_flag_default, grain=CLOCK_SUBCOMPONENT ) id_movnest4 = mpp_clock_id ('MN Part 4 Fill Intern Nest Halos', flags = clock_flag_default, grain=CLOCK_SUBCOMPONENT ) id_movnest5 = mpp_clock_id ('MN Part 5 Recalc Weights', flags = clock_flag_default, grain=CLOCK_SUBCOMPONENT ) + id_movnest5_1 = mpp_clock_id ('MN Part 5.1 read_parent', flags = clock_flag_default, grain=CLOCK_SUBCOMPONENT ) + id_movnest5_2 = mpp_clock_id ('MN Part 5.2 reset latlon', flags = clock_flag_default, grain=CLOCK_SUBCOMPONENT ) + id_movnest5_3 = mpp_clock_id ('MN Part 5.3 meta recalc', flags = clock_flag_default, grain=CLOCK_SUBCOMPONENT ) + id_movnest5_4 = mpp_clock_id ('MN Part 5.4 shift indx', flags = clock_flag_default, grain=CLOCK_SUBCOMPONENT ) + id_movnest6 = mpp_clock_id ('MN Part 6 EOSHIFT', flags = clock_flag_default, grain=CLOCK_SUBCOMPONENT ) id_movnest7_0 = mpp_clock_id ('MN Part 7.0 Recalc gridstruct', flags = clock_flag_default, grain=CLOCK_SUBCOMPONENT ) @@ -489,11 +495,14 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, type(grid_geometry), save :: parent_geo type(grid_geometry), save :: fp_super_tile_geo type(mn_surface_grids), save :: mn_static + real(kind=R_GRID), allocatable, save :: p_grid(:,:,:) + real(kind=R_GRID), allocatable, save :: p_grid_u(:,:,:) + real(kind=R_GRID), allocatable, save :: p_grid_v(:,:,:) type(grid_geometry) :: tile_geo, tile_geo_u, tile_geo_v - real(kind=R_GRID), allocatable :: p_grid(:,:,:), n_grid(:,:,:) - real(kind=R_GRID), allocatable :: p_grid_u(:,:,:), n_grid_u(:,:,:) - real(kind=R_GRID), allocatable :: p_grid_v(:,:,:), n_grid_v(:,:,:) + real(kind=R_GRID), allocatable :: n_grid(:,:,:) + real(kind=R_GRID), allocatable :: n_grid_u(:,:,:) + real(kind=R_GRID), allocatable :: n_grid_v(:,:,:) real, allocatable :: wt_h(:,:,:) ! TODO verify that these are deallocated real, allocatable :: wt_u(:,:,:) real, allocatable :: wt_v(:,:,:) @@ -913,17 +922,24 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, !!============================================================================ if (debug_log) print '("[INFO] WDR MV_NST5 run step 5 fv_moving_nest_main.F90 npe=",I0, " tile_geo%lats allocated:",L1)', this_pe, allocated(tile_geo%lats) if (debug_log) print '("[INFO] WDR MV_NST5 run step 5 fv_moving_nest_main.F90 npe=",I0, " parent_geo%lats allocated:",L1)', this_pe, allocated(parent_geo%lats) + if (use_timers) call mpp_clock_begin (id_movnest5_1) - ! parent_geo is only loaded first time; afterwards it is reused. - ! This is the coarse resolution data for the parent + ! parent_geo, p_grid, p_grid_u, and p_grid_v are only loaded first time; afterwards they are reused. + ! Because they are the coarse resolution grids (supergrid, a-grid, u stagger, v stagger) for the parent call mn_latlon_load_parent(Moving_nest(child_grid_num)%mn_flag%surface_dir, Atm, n, parent_tile, & delta_i_c, delta_j_c, child_grid_num, & parent_geo, tile_geo, tile_geo_u, tile_geo_v, fp_super_tile_geo, & p_grid, n_grid, p_grid_u, n_grid_u, p_grid_v, n_grid_v) + if (use_timers) call mpp_clock_end (id_movnest5_1) + if (use_timers) call mpp_clock_begin (id_movnest5_2) + ! tile_geo holds the center lat/lons for the entire nest (all PEs). 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) + !!============================================================================ !! Step 5.2 -- Fill the wt* variables for each stagger !!============================================================================ @@ -937,8 +953,11 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, 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) + if (use_timers) call mpp_clock_end (id_movnest5_3) endif + if (use_timers) call mpp_clock_begin (id_movnest5_4) + !!============================================================================ !! Step 5.3 -- Adjust the indices by the values of delta_i_c, delta_j_c !!============================================================================ @@ -950,6 +969,8 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, if (debug_sync) call mpp_sync(full_pelist) ! Used to make debugging easier. Can be removed. + if (use_timers) call mpp_clock_end (id_movnest5_4) + if (use_timers) call mpp_clock_end (id_movnest5) if (use_timers) call mpp_clock_begin (id_movnest6) diff --git a/moving_nest/fv_moving_nest_physics.F90 b/moving_nest/fv_moving_nest_physics.F90 index c2ff732fc..03a9dd38b 100644 --- a/moving_nest/fv_moving_nest_physics.F90 +++ b/moving_nest/fv_moving_nest_physics.F90 @@ -72,7 +72,6 @@ module fv_moving_nest_physics_mod use GFS_init, only: GFS_grid_populate use boundary_mod, only: update_coarse_grid, update_coarse_grid_mpp - use bounding_box_mod, only: bbox, bbox_get_C2F_index, fill_bbox, show_bbox use constants_mod, only: cp_air, rdgas, grav, rvgas, kappa, pstd_mks, hlv use field_manager_mod, only: MODEL_ATMOS use fms_io_mod, only: read_data, write_data, get_global_att_value, fms_io_init, fms_io_exit diff --git a/tools/fv_grid_tools.F90 b/tools/fv_grid_tools.F90 index 680d74fbf..ed3e56047 100644 --- a/tools/fv_grid_tools.F90 +++ b/tools/fv_grid_tools.F90 @@ -605,8 +605,8 @@ subroutine init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, ! Setup timing variables logical, save :: first_time = .true. - integer, save :: id_timer1, id_timer2, id_timer3, id_timer3a, id_timer4, id_timer5, id_timer6, id_timer7, id_timer8 - logical :: use_timer = .false. ! Set to True for detailed performance profiling + integer, save :: id_timer1, id_timer2, id_timer3, id_timer3a, id_timer3b, id_timer4, id_timer5, id_timer6, id_timer7, id_timer8 + logical :: use_timer = .True. ! Set to True for detailed performance profiling logical :: debug_log = .false. integer :: this_pe @@ -617,7 +617,8 @@ subroutine init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, id_timer1 = mpp_clock_id ('init_grid Step 1', flags = clock_flag_default, grain=CLOCK_ROUTINE ) id_timer2 = mpp_clock_id ('init_grid Step 2', flags = clock_flag_default, grain=CLOCK_ROUTINE ) id_timer3 = mpp_clock_id ('init_grid Step 3', flags = clock_flag_default, grain=CLOCK_ROUTINE ) - id_timer3a = mpp_clock_id ('init_grid Step 3a', flags = clock_flag_default, grain=CLOCK_ROUTINE ) + id_timer3a = mpp_clock_id ('init_grid Step 3a read_grid', flags = clock_flag_default, grain=CLOCK_ROUTINE ) + id_timer3b = mpp_clock_id ('init_grid Step 3b setup_aligned_nest', flags = clock_flag_default, grain=CLOCK_ROUTINE ) id_timer4 = mpp_clock_id ('init_grid Step 4', flags = clock_flag_default, grain=CLOCK_ROUTINE ) id_timer5 = mpp_clock_id ('init_grid Step 5', flags = clock_flag_default, grain=CLOCK_ROUTINE ) id_timer6 = mpp_clock_id ('init_grid Step 6', flags = clock_flag_default, grain=CLOCK_ROUTINE ) @@ -741,15 +742,17 @@ subroutine init_grid(Atm, grid_name, grid_file, npx, npy, npz, ndims, nregions, if (Atm%neststruct%nested) then !Read grid if it exists - if (use_timer) call mpp_clock_begin (id_timer3a) if (Atm%flagstruct%grid_type < 0) then + if (use_timer) call mpp_clock_begin (id_timer3a) !Note that read_grid only reads in grid corners. Will still need to compute all other grid metrics. !NOTE: cannot currently read in mosaic for both coarse and nested grids simultaneously call read_grid(Atm, grid_file, ndims, 1, ng) + if (use_timer) call mpp_clock_end (id_timer3a) endif ! still need to set up weights + if (use_timer) call mpp_clock_begin (id_timer3b) call setup_aligned_nest(Atm) - if (use_timer) call mpp_clock_end (id_timer3a) + if (use_timer) call mpp_clock_end (id_timer3b) else if(trim(grid_file) .NE. 'Inline' .or. Atm%flagstruct%grid_type < 0) then From 0cd5847e4de0cb48b0ac27d450e7db4b14df62c5 Mon Sep 17 00:00:00 2001 From: "Bin.Liu" Date: Fri, 20 May 2022 19:23:38 -0500 Subject: [PATCH 03/29] Update atmos_model_nml in driver/fvGFS/atmosphere.F90 so that it is consistent with that in FV3/atmos_model.F90. --- driver/fvGFS/atmosphere.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/driver/fvGFS/atmosphere.F90 b/driver/fvGFS/atmosphere.F90 index f2f68eb6e..86c0fda42 100644 --- a/driver/fvGFS/atmosphere.F90 +++ b/driver/fvGFS/atmosphere.F90 @@ -324,10 +324,10 @@ subroutine atmosphere_init (Time_init, Time, Time_step, Grid_box, area) logical :: dycore_only = .false. logical :: debug = .false. logical :: sync = .false. - integer, parameter :: maxhr = 4096 - real, dimension(maxhr) :: fdiag = 0. - real :: fhmax=384.0, fhmaxhf=120.0, fhout=3.0, fhouthf=1.0,avg_max_length=3600. - namelist /atmos_model_nml/ blocksize, chksum_debug, dycore_only, debug, sync, fdiag, fhmax, fhmaxhf, fhout, fhouthf, ccpp_suite, avg_max_length + real :: avg_max_length=3600. + logical :: ignore_rst_cksum = .false. + namelist /atmos_model_nml/ blocksize, chksum_debug, dycore_only, debug, sync, ccpp_suite, avg_max_length, & + ignore_rst_cksum ! *DH 20210326 !For regional From 205d7752a5e1785ae171716dda04c5374b42115a Mon Sep 17 00:00:00 2001 From: William Ramstrom Date: Fri, 29 Jul 2022 18:31:57 +0000 Subject: [PATCH 04/29] Removed reference to unused variable parent_x. --- moving_nest/fv_moving_nest.F90 | 7 ------- 1 file changed, 7 deletions(-) diff --git a/moving_nest/fv_moving_nest.F90 b/moving_nest/fv_moving_nest.F90 index 1e1f2a8be..252af6291 100644 --- a/moving_nest/fv_moving_nest.F90 +++ b/moving_nest/fv_moving_nest.F90 @@ -681,7 +681,6 @@ subroutine mn_latlon_load_parent(surface_dir, Atm, n, parent_tile, delta_i_c, de integer :: x, y, fp_i, fp_j integer :: position, position_u, position_v integer :: x_refine, y_refine - integer :: nest_x, nest_y, parent_x, parent_y integer :: this_pe logical, save :: first_time = .True. @@ -852,12 +851,6 @@ subroutine mn_latlon_load_parent(surface_dir, Atm, n, parent_tile, delta_i_c, de if (use_timers) call mpp_clock_end (id_load4) if (use_timers) call mpp_clock_begin (id_load5) - - if (parent_x .eq. -999) then - print '("[ERROR] WDR mn_latlon_load_parent on npe=",I0," parent and nest grids are not aligned!")', this_pe - call mpp_error(FATAL, "mn_latlon_load_parent parent and nest grids are not aligned.") - endif - ! These grids are small (nest size), and change each time nest moves. call assign_n_grids(tile_geo, n_grid, position) call assign_n_grids(tile_geo_u, n_grid_u, position_u) From af29d30627c49555199096c058914af1b83dacb3 Mon Sep 17 00:00:00 2001 From: Biju Thomas Date: Wed, 12 Oct 2022 13:32:36 +0000 Subject: [PATCH 05/29] Adding upoff as a namelist parameter --- model/fv_control.F90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/model/fv_control.F90 b/model/fv_control.F90 index 6951a15a5..f03edb69d 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -408,7 +408,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split, integer, pointer :: nrows_blend logical, pointer :: regional_bcs_from_gsi logical, pointer :: write_restart_with_bcs - integer, pointer :: parent_tile, refinement, nestbctype, nestupdate, nsponge, ioffset, joffset + integer, pointer :: parent_tile, refinement, nestbctype, nestupdate, upoff, nsponge, ioffset, joffset real, pointer :: s_weight, update_blend character(len=16), pointer :: restart_resolution @@ -971,6 +971,7 @@ subroutine set_namelist_pointers(Atm) refinement => Atm%neststruct%refinement nestbctype => Atm%neststruct%nestbctype nestupdate => Atm%neststruct%nestupdate + upoff => Atm%neststruct%upoff nsponge => Atm%neststruct%nsponge s_weight => Atm%neststruct%s_weight ioffset => Atm%neststruct%ioffset @@ -1071,7 +1072,7 @@ subroutine read_namelist_fv_core_nml(Atm) deglon_start, deglon_stop, deglat_start, deglat_stop, & phys_hydrostatic, use_hydro_pressure, make_hybrid_z, old_divg_damp, add_noise, butterfly_effect, & molecular_diffusion, dz_min, psm_bc, nested, twowaynest, nudge_qv, & - nestbctype, nestupdate, nsponge, s_weight, & + nestbctype, nestupdate, upoff, nsponge, s_weight, & check_negative, nudge_ic, halo_update_type, gfs_phil, agrid_vel_rst, & do_uni_zfull, adj_mass_vmr, fac_n_spl, fhouri, update_blend, regional, bc_update_interval, & regional_bcs_from_gsi, write_restart_with_bcs, nrows_blend, & From 9b56598dc1f3a50a42eca32268b2e7cbf1967936 Mon Sep 17 00:00:00 2001 From: Bin Liu Date: Fri, 28 Oct 2022 17:14:57 +0000 Subject: [PATCH 06/29] Update to output timestr in yyyymmdd.hhmmss for the internal tracker output (fort.602, the partial atcf track file). --- moving_nest/fv_tracker.F90 | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/moving_nest/fv_tracker.F90 b/moving_nest/fv_tracker.F90 index f00c9da57..40780864d 100644 --- a/moving_nest/fv_tracker.F90 +++ b/moving_nest/fv_tracker.F90 @@ -30,7 +30,7 @@ module fv_tracker_mod use constants_mod, only: pi=>pi_8, rad_to_deg, deg_to_rad use time_manager_mod, only: time_type, get_time, set_time, operator(+), & - operator(-), operator(/), time_type_to_real + operator(-), operator(/), time_type_to_real, date_to_string use mpp_mod, only: mpp_error, stdout, FATAL, WARNING, NOTE, & mpp_root_pe, mpp_npes, mpp_pe, mpp_chksum, & mpp_get_current_pelist, & @@ -744,24 +744,23 @@ subroutine output_partial_atcfunix(Atm,Time, & integer, intent(in) :: ids,ide,jds,jde,kds,kde integer, intent(in) :: ims,ime,jms,jme,kms,kme integer, intent(in) :: its,ite,jts,jte,kts,kte - integer :: days, seconds - real :: sec + character*15 timestr character*255 message - call get_time(fv_time, seconds, days) - sec=seconds -313 format(F11.2,", ", & + ! timestr in the format of yyyymmdd.hhmmss + timestr=date_to_string(fv_time) +313 format(A15,", ", & "W10 = ",F7.3," kn, PMIN = ",F8.3," mbar, ", & "LAT = ",F6.3,A1,", LON = ",F7.3,A1,", ", & "RMW = ",F7.3," nmi") if (Tracker(n)%tracker_fixlon .gt. 180.0) then - write(Moving_nest(n)%mn_flag%outatcf_lun+Atm%grid_number,313) sec, & + write(Moving_nest(n)%mn_flag%outatcf_lun+Atm%grid_number,313) timestr, & Tracker(n)%tracker_vmax*mps2kn,Tracker(n)%tracker_pmin/100., & abs(Tracker(n)%tracker_fixlat),get_lat_ns(Tracker(n)%tracker_fixlat), & abs(Tracker(n)%tracker_fixlon-360.0),get_lon_ew(Tracker(n)%tracker_fixlon-360.0), & Tracker(n)%tracker_rmw*km2nmi else - write(Moving_nest(n)%mn_flag%outatcf_lun+Atm%grid_number,313) sec, & + write(Moving_nest(n)%mn_flag%outatcf_lun+Atm%grid_number,313) timestr, & Tracker(n)%tracker_vmax*mps2kn,Tracker(n)%tracker_pmin/100., & abs(Tracker(n)%tracker_fixlat),get_lat_ns(Tracker(n)%tracker_fixlat), & abs(Tracker(n)%tracker_fixlon),get_lon_ew(Tracker(n)%tracker_fixlon), & From 59a74c9cccc9864ff9efbcffdc5df40af9914cf6 Mon Sep 17 00:00:00 2001 From: Bin Liu Date: Thu, 23 Feb 2023 20:17:25 +0000 Subject: [PATCH 07/29] Remove "-nowarn" build option in cmake/compiler_flags_Intel_Fortran.cmake (from @BijuThomas-NOAA). --- cmake/compiler_flags_Intel_Fortran.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/compiler_flags_Intel_Fortran.cmake b/cmake/compiler_flags_Intel_Fortran.cmake index 743afce35..0cfc70397 100644 --- a/cmake/compiler_flags_Intel_Fortran.cmake +++ b/cmake/compiler_flags_Intel_Fortran.cmake @@ -4,7 +4,7 @@ set(R8_flags "-real-size 64") # Fortran flags for 64BIT precision set(R8_flags "${R8_flags} -no-prec-div -no-prec-sqrt") # Intel Fortran -set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -g -traceback -fpp -fno-alias -auto -safe-cray-ptr -ftz -assume byterecl -nowarn -sox -align array64byte -qno-opt-dynamic-align ${${kind}_flags}") +set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -g -traceback -fpp -fno-alias -auto -safe-cray-ptr -ftz -assume byterecl -sox -align array64byte -qno-opt-dynamic-align ${${kind}_flags}") set(CMAKE_Fortran_FLAGS_REPRO "-O2 -debug minimal -fp-model consistent -qoverride-limits") From cc648261d6c9846bf1ccb3b88eaa1b145b7d1190 Mon Sep 17 00:00:00 2001 From: Dusan Jovic <48258889+DusanJovic-NOAA@users.noreply.github.com> Date: Wed, 25 Jan 2023 09:24:32 -0500 Subject: [PATCH 08/29] Use 32bit value for 'missing_value' and '_FillValue' attributes (#236) --- driver/fvGFS/fv_nggps_diag.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/driver/fvGFS/fv_nggps_diag.F90 b/driver/fvGFS/fv_nggps_diag.F90 index d2b64a5d9..3d05ecfe5 100644 --- a/driver/fvGFS/fv_nggps_diag.F90 +++ b/driver/fvGFS/fv_nggps_diag.F90 @@ -1727,12 +1727,12 @@ subroutine add_field_to_bundle(var_name,long_name,units,cell_methods, axes,dyn_g call ESMF_AttributeAdd(field, convention="NetCDF", purpose="FV3", & attrList=(/"missing_value"/), rc=rc) call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3", & - name='missing_value',value=missing_value,rc=rc) + name='missing_value',value=real(missing_value,kind=4),rc=rc) call ESMF_AttributeAdd(field, convention="NetCDF", purpose="FV3", & attrList=(/"_FillValue"/), rc=rc) call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3", & - name='_FillValue',value=missing_value,rc=rc) + name='_FillValue',value=real(missing_value,kind=4),rc=rc) call ESMF_AttributeAdd(field, convention="NetCDF", purpose="FV3", & attrList=(/"cell_methods"/), rc=rc) From 94bf82f5456ad6a56ee2e016a6387996558935ba Mon Sep 17 00:00:00 2001 From: Rusty Benson Date: Tue, 27 Feb 2024 12:16:12 -0500 Subject: [PATCH 09/29] give istatus a value on exit from hailcast_init --- tools/module_diag_hailcast.F90 | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tools/module_diag_hailcast.F90 b/tools/module_diag_hailcast.F90 index 2a8e37a33..0a9e78525 100644 --- a/tools/module_diag_hailcast.F90 +++ b/tools/module_diag_hailcast.F90 @@ -84,6 +84,11 @@ SUBROUTINE hailcast_init(file_name, axes, Time, isco,ieco,jsco,jeco,& !write(unit, nml=fv_diagnostics_nml) !!end hailcast nml + + ! need to set a default value for istatus because of it's intent(out) status + ! this value is not checked on return + istatus = 0 + if (mpp_pe() == mpp_root_pe()) then print*, 'do_hailcast = ', do_hailcast end if From f7ff2afd4cde8af81f79603f1ccb667dae8b480a Mon Sep 17 00:00:00 2001 From: Rusty Benson Date: Tue, 27 Feb 2024 13:58:00 -0500 Subject: [PATCH 10/29] fix issues with declared out variables in test_cases.F90 --- tools/test_cases.F90 | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/tools/test_cases.F90 b/tools/test_cases.F90 index 420958149..f49f1cb25 100644 --- a/tools/test_cases.F90 +++ b/tools/test_cases.F90 @@ -3192,7 +3192,7 @@ subroutine init_case(u,v,w,pt,delp,q,phis, ps,pe,peln,pk,pkz, uc,vc, ua,va, ak, ! Iterate then interpolate to get balanced pt & pk on the sphere ! Adjusting ptop call SuperK_u(npz, zs1, uz1, dudz) - call balanced_K(npz, is, ie, js, je, ng, pe1(npz+1), ze1, ts1, qs1, uz1, dudz, pe, pk, pt, & + call balanced_K(npz, is, ie, js, je, ng, pe1(npz+1), ze1, ts1, qs1, uz1, dudz, pe, pt, & delz, zvir, ptop, ak, bk, agrid) do j=js,je do i=is,ie @@ -5464,7 +5464,7 @@ subroutine SuperK_Sounding(km, pe, p00, ze, pt, qz) end subroutine SuperK_Sounding - subroutine balanced_K(km, is, ie, js, je, ng, ps0, ze1, ts1, qs1, uz1, dudz, pe, pk, pt, & + subroutine balanced_K(km, is, ie, js, je, ng, ps0, ze1, ts1, qs1, uz1, dudz, pe, pt, & delz, zvir, ptop, ak, bk, agrid) integer, intent(in):: is, ie, js, je, ng, km real, intent(in), dimension(km ):: ts1, qs1, uz1, dudz @@ -5475,7 +5475,6 @@ subroutine balanced_K(km, is, ie, js, je, ng, ps0, ze1, ts1, qs1, uz1, dudz, pe, real, intent(inout), dimension(km+1):: ak, bk real, intent(inout), dimension(is:ie,js:je,km):: pt real, intent(inout), dimension(is:,js:,1:) :: delz - real, intent(out), dimension(is:ie,js:je,km+1):: pk ! pt is FV's cp*thelta_v real, intent(inout), dimension(is-1:ie+1,km+1,js-1:je+1):: pe ! Local @@ -5683,6 +5682,8 @@ subroutine SuperCell_Sounding(km, ps, pk1, tp, qp) #ifdef GFS_PHYS call mpp_error(FATAL, 'SuperCell sounding cannot perform with GFS Physics.') + tp=0. + qb=0. #else From 978f1def59b474da10d2376648caf4bda22b5810 Mon Sep 17 00:00:00 2001 From: Rusty Benson Date: Wed, 28 Feb 2024 14:05:20 -0500 Subject: [PATCH 11/29] fix out of range default integers, eol spaces, and remove delz argument from set_regional_BCs in fv_regional_bc.F90 --- model/fv_dynamics.F90 | 2 +- model/fv_regional_bc.F90 | 36 ++++++++++++++++++------------------ 2 files changed, 19 insertions(+), 19 deletions(-) diff --git a/model/fv_dynamics.F90 b/model/fv_dynamics.F90 index 31aaa62a1..d4881c23d 100644 --- a/model/fv_dynamics.F90 +++ b/model/fv_dynamics.F90 @@ -378,7 +378,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, reg_bc_update_time=current_time_in_seconds call set_regional_BCs & !<-- Insert values into the boundary region valid for the start of this large timestep. - (delp,delz,w,pt & + (delp,w,pt & #ifdef USE_COND ,q_con & #endif diff --git a/model/fv_regional_bc.F90 b/model/fv_regional_bc.F90 index 8b672e1fa..ba098334d 100644 --- a/model/fv_regional_bc.F90 +++ b/model/fv_regional_bc.F90 @@ -93,6 +93,7 @@ module fv_regional_mod integer,parameter :: bc_time_interval=3 & ,nhalo_data =4 & ,nhalo_model=3 + integer, public, parameter :: int_init_default = -9999999 ! integer, public, parameter :: H_STAGGER = 1 integer, public, parameter :: U_STAGGER = 2 @@ -471,7 +472,7 @@ subroutine setup_regional_BC(Atm & else nrows_blend=nrows_blend_in_data !<-- # of blending rows in the BC files. endif - + IF ( north_bc .or. south_bc ) THEN IF ( nrows_blend_user > jed - nhalo_model - (jsd + nhalo_model) + 1 ) THEN call mpp_error(FATAL,'Number of blending rows is greater than the north-south tile size!') @@ -4076,7 +4077,7 @@ end subroutine remap_dwinds_regional_bc !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !--------------------------------------------------------------------- - subroutine set_regional_BCs(delp,delz,w,pt & + subroutine set_regional_BCs(delp,w,pt & #ifdef USE_COND ,q_con & #endif @@ -4085,7 +4086,7 @@ subroutine set_regional_BCs(delp,delz,w,pt & #endif ,q & ,u,v,uc,vc & - ,bd, nlayers & + ,bd, nlayers & ,fcst_time ) ! !--------------------------------------------------------------------- @@ -4117,7 +4118,6 @@ subroutine set_regional_BCs(delp,delz,w,pt & ,pt ! real,dimension(bd%isd:,bd%jsd:,1:),intent(out) :: w - real,dimension(bd%is:,bd%js:,1:),intent(out) :: delz #ifdef USE_COND real,dimension(bd%isd:,bd%jsd:,1:),intent(out) :: q_con #endif @@ -4404,7 +4404,7 @@ subroutine regional_boundary_update(array & ! integer,intent(in) :: is,ie,js,je & !<-- Compute limits ,isd,ied,jsd,jed & !<-- Memory limits - ,it !<-- Acoustic step + ,it !<-- Acoustic step ! integer,intent(in),optional :: index4 !<-- Index for the 4-D tracer array. ! @@ -4494,7 +4494,7 @@ subroutine regional_boundary_update(array & endif j1_blend=js j2_blend=js+nrows_blend_user-1 - i_bc=-9e9 + i_bc=int_init_default j_bc=j2 ! endif @@ -4544,7 +4544,7 @@ subroutine regional_boundary_update(array & j2_blend=je+1 endif j1_blend=j2_blend-nrows_blend_user+1 - i_bc=-9e9 + i_bc=int_init_default j_bc=j1 ! endif @@ -4601,7 +4601,7 @@ subroutine regional_boundary_update(array & j2_blend=j2_blend+1 endif i_bc=i2 - j_bc=-9e9 + j_bc=int_init_default ! endif endif @@ -4660,7 +4660,7 @@ subroutine regional_boundary_update(array & j2_blend=j2_blend+1 endif i_bc=i1 - j_bc=-9e9 + j_bc=int_init_default ! endif endif @@ -6892,10 +6892,10 @@ subroutine get_data_source(data_source_fv3gfs,regional) if (.not. lstatus) then if (mpp_pe() == 0) write(0,*) 'INPUT source not found ',lstatus,' set source=No Source Attribute' source='No Source Attribute' - call mpp_error(FATAL,'fv_regional_bc::get_data_source - input source not & - found in file gfs_data.nc. The accepted & + call mpp_error(FATAL,'fv_regional_bc::get_data_source - input source not & + found in file gfs_data.nc. The accepted & FV3 sources are "FV3GFS GAUSSIAN NEMSIO FILE", & - "FV3GFS GAUSSIAN NETCDF FILE" or "FV3GFS GRIB2 FILE".') + "FV3GFS GAUSSIAN NETCDF FILE" or "FV3GFS GRIB2 FILE".') endif call mpp_error(NOTE, 'INPUT gfs_data source string: '//trim(source)) @@ -6925,7 +6925,7 @@ subroutine get_lbc_source(lbc_source_fv3gfs,regional) character (len=80) :: source logical :: lstatus = .false. type(FmsNetcdfFile_t) :: Gfs_data - integer, allocatable, dimension(:) :: pes !< Array of the pes in the current pelist + integer, allocatable, dimension(:) :: pes !< Array of the pes in the current pelist ! ! Use the fms call here so we can actually get the return code value. ! The term 'source' is specified by 'chgres_cube' @@ -6934,7 +6934,7 @@ subroutine get_lbc_source(lbc_source_fv3gfs,regional) allocate(pes(mpp_npes())) call mpp_get_current_pelist(pes) - if (open_file(Gfs_data , 'INPUT/gfs_bndy.tile7.000.nc', "read", pelist=pes)) then + if (open_file(Gfs_data , 'INPUT/gfs_bndy.tile7.000.nc', "read", pelist=pes)) then lstatus = global_att_exists(Gfs_data, "source") if(lstatus) call get_global_attribute(Gfs_data, "source", source) call close_file(Gfs_data) @@ -6942,13 +6942,13 @@ subroutine get_lbc_source(lbc_source_fv3gfs,regional) deallocate(pes) if (.not. lstatus) then - if (mpp_pe() == 0) write(0,*) 'INPUT source not found ',lstatus,' set source=No Source Attribute' + if (mpp_pe() == 0) write(0,*) 'INPUT source not found ',lstatus,' set source=No Source Attribute' source='No Source Attribute' - call mpp_error(FATAL,'fv_regional_bc::get_lbc_source - input source not & + call mpp_error(FATAL,'fv_regional_bc::get_lbc_source - input source not & found in file & - gfs_bndy.tile7.000.nc. The accepted & + gfs_bndy.tile7.000.nc. The accepted & FV3 sources are "FV3GFS GAUSSIAN NEMSIO FILE", & - "FV3GFS GAUSSIAN NETCDF FILE" or "FV3GFS GRIB2 FILE".') + "FV3GFS GAUSSIAN NETCDF FILE" or "FV3GFS GRIB2 FILE".') endif call mpp_error(NOTE, 'INPUT gfs_bndy source string: '//trim(source)) From 2cc9f5d76aa5ccaf6c4821df9c39b5ae4c78aca5 Mon Sep 17 00:00:00 2001 From: Rusty Benson Date: Wed, 28 Feb 2024 15:04:45 -0500 Subject: [PATCH 12/29] removed -nowarn from Intel compiler flags --- cmake/compiler_flags_Intel_Fortran.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/compiler_flags_Intel_Fortran.cmake b/cmake/compiler_flags_Intel_Fortran.cmake index 743afce35..0cfc70397 100644 --- a/cmake/compiler_flags_Intel_Fortran.cmake +++ b/cmake/compiler_flags_Intel_Fortran.cmake @@ -4,7 +4,7 @@ set(R8_flags "-real-size 64") # Fortran flags for 64BIT precision set(R8_flags "${R8_flags} -no-prec-div -no-prec-sqrt") # Intel Fortran -set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -g -traceback -fpp -fno-alias -auto -safe-cray-ptr -ftz -assume byterecl -nowarn -sox -align array64byte -qno-opt-dynamic-align ${${kind}_flags}") +set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -g -traceback -fpp -fno-alias -auto -safe-cray-ptr -ftz -assume byterecl -sox -align array64byte -qno-opt-dynamic-align ${${kind}_flags}") set(CMAKE_Fortran_FLAGS_REPRO "-O2 -debug minimal -fp-model consistent -qoverride-limits") From 4bfb82d283e8bc5d025d3988a83183be9f24ee35 Mon Sep 17 00:00:00 2001 From: Rusty Benson Date: Wed, 28 Feb 2024 15:07:09 -0500 Subject: [PATCH 13/29] typo qb->qp --- tools/test_cases.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/test_cases.F90 b/tools/test_cases.F90 index f49f1cb25..49f86589e 100644 --- a/tools/test_cases.F90 +++ b/tools/test_cases.F90 @@ -5683,7 +5683,7 @@ subroutine SuperCell_Sounding(km, ps, pk1, tp, qp) call mpp_error(FATAL, 'SuperCell sounding cannot perform with GFS Physics.') tp=0. - qb=0. + qp=0. #else From 6bcd3bb249399004649f865b8fa73b6431f369e3 Mon Sep 17 00:00:00 2001 From: Rusty Benson Date: Fri, 29 Mar 2024 10:05:23 -0400 Subject: [PATCH 14/29] added suggested ifdef to remove goto warning as per operations --- tools/fv_nudge.F90 | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/tools/fv_nudge.F90 b/tools/fv_nudge.F90 index 787edefa9..f79f5f919 100644 --- a/tools/fv_nudge.F90 +++ b/tools/fv_nudge.F90 @@ -2273,6 +2273,10 @@ subroutine breed_slp_inline(nstep, dt, npz, ak, bk, phis, pe, pk, peln, pkz, del real :: kappax(is:ie,js:je,npz) #endif +#if defined (BYPASS_BREED_SLP_INLINE) + peln = 0.0 ! to silence compiler warning. A dummy argument with an explicit INTENT(OUT) declaration is not given an explicit value. + call mpp_error(fatal, "breed_slp_inline routine has been disabled") +#else if ( forecast_mode ) return agrid => gridstruct%agrid_64 @@ -2715,6 +2719,7 @@ subroutine breed_slp_inline(nstep, dt, npz, ak, bk, phis, pe, pk, peln, pkz, del nullify(agrid) nullify(area) +#endif end subroutine breed_slp_inline From 6fcbe1357dff6b5bdc4abd67c342b56ab2de935c Mon Sep 17 00:00:00 2001 From: BijuThomas-NOAA Date: Fri, 5 Apr 2024 13:58:26 +0000 Subject: [PATCH 15/29] Remove the nowarn flag from the build --- cmake/compiler_flags_Intel_Fortran.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/compiler_flags_Intel_Fortran.cmake b/cmake/compiler_flags_Intel_Fortran.cmake index 743afce35..0cfc70397 100644 --- a/cmake/compiler_flags_Intel_Fortran.cmake +++ b/cmake/compiler_flags_Intel_Fortran.cmake @@ -4,7 +4,7 @@ set(R8_flags "-real-size 64") # Fortran flags for 64BIT precision set(R8_flags "${R8_flags} -no-prec-div -no-prec-sqrt") # Intel Fortran -set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -g -traceback -fpp -fno-alias -auto -safe-cray-ptr -ftz -assume byterecl -nowarn -sox -align array64byte -qno-opt-dynamic-align ${${kind}_flags}") +set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -g -traceback -fpp -fno-alias -auto -safe-cray-ptr -ftz -assume byterecl -sox -align array64byte -qno-opt-dynamic-align ${${kind}_flags}") set(CMAKE_Fortran_FLAGS_REPRO "-O2 -debug minimal -fp-model consistent -qoverride-limits") From 4162fad68d84db507bf25d4502673bc608c86600 Mon Sep 17 00:00:00 2001 From: Bin Liu Date: Mon, 8 Apr 2024 20:45:48 +0000 Subject: [PATCH 16/29] Correct the jc index when print out the nest grid corner locations in tools/fv_grid_tools.F90. --- tools/fv_grid_tools.F90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/tools/fv_grid_tools.F90 b/tools/fv_grid_tools.F90 index d74000f64..67b33ccf2 100644 --- a/tools/fv_grid_tools.F90 +++ b/tools/fv_grid_tools.F90 @@ -2655,23 +2655,23 @@ subroutine setup_aligned_nest(Atm) !Nesting position information !BUG multiply by 180 not 90.... write(*,*) 'NESTED GRID ', Atm%grid_number - ic = p_ind(1,1,1) ; jc = p_ind(1,1,1) + ic = p_ind(1,1,1) ; jc = p_ind(1,1,2) write(*,'(A, 2I5, 4F10.4)') 'SW CORNER: ', ic, jc, grid_global(1,1,:,1)*90./pi - ic = p_ind(1,npy,1) ; jc = p_ind(1,npy,1) + ic = p_ind(1,npy,1) ; jc = p_ind(1,npy,2) write(*,'(A, 2I5, 4F10.4)') 'NW CORNER: ', ic, jc, grid_global(1,npy,:,1)*90./pi - ic = p_ind(npx,npy,1) ; jc = p_ind(npx,npy,1) + ic = p_ind(npx,npy,1) ; jc = p_ind(npx,npy,2) write(*,'(A, 2I5, 4F10.4)') 'NE CORNER: ', ic, jc, grid_global(npx,npy,:,1)*90./pi - ic = p_ind(npx,1,1) ; jc = p_ind(npx,1,1) + ic = p_ind(npx,1,1) ; jc = p_ind(npx,1,2) write(*,'(A, 2I5, 4F10.4)') 'SE CORNER: ', ic, jc, grid_global(npx,1,:,1)*90./pi else write(*,*) 'PARENT GRID ', Atm%parent_grid%grid_number, Atm%parent_grid%global_tile - ic = p_ind(1,1,1) ; jc = p_ind(1,1,1) + ic = p_ind(1,1,1) ; jc = p_ind(1,1,2) write(*,'(A, 2I5, 4F10.4)') 'SW CORNER: ', ic, jc, Atm%parent_grid%grid_global(ic,jc,:,parent_tile)*90./pi - ic = p_ind(1,npy,1) ; jc = p_ind(1,npy,1) + ic = p_ind(1,npy,1) ; jc = p_ind(1,npy,2) write(*,'(A, 2I5, 4F10.4)') 'NW CORNER: ', ic, jc, Atm%parent_grid%grid_global(ic,jc,:,parent_tile)*90./pi - ic = p_ind(npx,npy,1) ; jc = p_ind(npx,npy,1) + ic = p_ind(npx,npy,1) ; jc = p_ind(npx,npy,2) write(*,'(A, 2I5, 4F10.4)') 'NE CORNER: ', ic, jc, Atm%parent_grid%grid_global(ic,jc,:,parent_tile)*90./pi - ic = p_ind(npx,1,1) ; jc = p_ind(npx,1,1) + ic = p_ind(npx,1,1) ; jc = p_ind(npx,1,2) write(*,'(A, 2I5, 4F10.4)') 'SE CORNER: ', ic, jc, Atm%parent_grid%grid_global(ic,jc,:,parent_tile)*90./pi endif end if From 3bdb6a1cfc2894814c5a9473963f6759eff3ecaf Mon Sep 17 00:00:00 2001 From: XuLu-NOAA Date: Fri, 12 Jul 2024 18:19:28 +0000 Subject: [PATCH 17/29] Add IAU related functions for regional HAFS DA. --- tools/fv_iau_mod.F90 | 97 ++++++++++++++++++++++++++++++++------------ 1 file changed, 70 insertions(+), 27 deletions(-) diff --git a/tools/fv_iau_mod.F90 b/tools/fv_iau_mod.F90 index e4ca56500..9db598e19 100644 --- a/tools/fv_iau_mod.F90 +++ b/tools/fv_iau_mod.F90 @@ -188,9 +188,15 @@ subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm) if( file_exists(fname) ) then call open_ncfile( fname, ncid ) ! open the file - call get_ncdim1( ncid, 'lon', im) - call get_ncdim1( ncid, 'lat', jm) - call get_ncdim1( ncid, 'lev', km) + if (IPD_Control%iau_regional) then + call get_ncdim1( ncid, 'nx', im) + call get_ncdim1( ncid, 'ny', jm) + call get_ncdim1( ncid, 'nz', km) + else + call get_ncdim1( ncid, 'lon', im) + call get_ncdim1( ncid, 'lat', jm) + call get_ncdim1( ncid, 'lev', km) + end if if (km.ne.npz) then if (is_master()) print *, 'km = ', km @@ -200,20 +206,22 @@ subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm) if(is_master()) write(*,*) fname, ' DA increment dimensions:', im,jm,km - allocate ( lon(im) ) - allocate ( lat(jm) ) + if (.not.IPD_Control%iau_regional) then + allocate ( lon(im) ) + allocate ( lat(jm) ) - call _GET_VAR1 (ncid, 'lon', im, lon ) - call _GET_VAR1 (ncid, 'lat', jm, lat ) - call close_ncfile(ncid) + call _GET_VAR1 (ncid, 'lon', im, lon ) + call _GET_VAR1 (ncid, 'lat', jm, lat ) + call close_ncfile(ncid) - ! Convert to radians - do i=1,im - lon(i) = lon(i) * deg2rad - enddo - do j=1,jm - lat(j) = lat(j) * deg2rad - enddo + ! Convert to radians + do i=1,im + lon(i) = lon(i) * deg2rad + enddo + do j=1,jm + lat(j) = lat(j) * deg2rad + enddo + end if else call mpp_error(FATAL,'==> Error in IAU_initialize: Expected file '& @@ -232,10 +240,12 @@ subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm) agrid(is-1+i,js-1+j,2)=Init_parm%xlat(i,j) enddo enddo - call remap_coef( is, ie, js, je, is, ie, js, je, & + if (.not.IPD_Control%iau_regional) then + call remap_coef( is, ie, js, je, is, ie, js, je, & im, jm, lon, lat, id1, id2, jdc, s2c, & agrid) - deallocate ( lon, lat,agrid ) + deallocate ( lon, lat,agrid ) + end if allocate(IAU_Data%ua_inc(is:ie, js:je, km)) @@ -253,6 +263,7 @@ subroutine IAU_initialize (IPD_Control, IAU_Data,Init_parm) allocate (iau_state%inc1%tracer_inc(is:ie, js:je, km,ntracers)) iau_state%hr1=IPD_Control%iaufhrs(1) iau_state%wt = 1.0 ! IAU increment filter weights (default 1.0) + iau_state%wt=iau_state%wt*IPD_Control%iau_inc_scale !Increase the weight iau_state%wt_normfact = 1.0 if (IPD_Control%iau_filter_increments) then ! compute increment filter weights, sum to obtain normalization factor @@ -346,7 +357,7 @@ subroutine getiauforcing(IPD_Control,IAU_Data) IAU_Data%in_interval=.false. else if (IPD_Control%iau_filter_increments) call setiauforcing(IPD_Control,IAU_Data,iau_state%wt) - if (is_master()) print *,'apply iau forcing t1,t,t2,filter wt=',t1,IPD_Control%fhour,t2,iau_state%wt/iau_state%wt_normfact + if (is_master()) print *,'apply iau forcing t1,t,t2,filter wt=',t1,IPD_Control%fhour,t2,iau_state%wt/iau_state%wt_normfact,iau_state%wt IAU_Data%in_interval=.true. endif return @@ -474,16 +485,32 @@ subroutine read_iau_forcing(IPD_Control,increments,fname) enddo enddo - allocate ( wk3(1:im,jbeg:jend, 1:km) ) + if (IPD_Control%iau_regional) then + allocate ( wk3(1:im,js:je, 1:km) ) + else + allocate ( wk3(1:im,jbeg:jend, 1:km) ) + endif + ! read in 1 time level - call interp_inc('T_inc',increments%temp_inc(:,:,:),jbeg,jend) - call interp_inc('delp_inc',increments%delp_inc(:,:,:),jbeg,jend) - call interp_inc('delz_inc',increments%delz_inc(:,:,:),jbeg,jend) - call interp_inc('u_inc',increments%ua_inc(:,:,:),jbeg,jend) ! can these be treated as scalars? - call interp_inc('v_inc',increments%va_inc(:,:,:),jbeg,jend) - do l=1,ntracers - call interp_inc(trim(tracer_names(l))//'_inc',increments%tracer_inc(:,:,:,l),jbeg,jend) - enddo + if (IPD_Control%iau_regional) then + call interp_inc2('T_inc',increments%temp_inc) + call interp_inc2('delp_inc',increments%delp_inc) + call interp_inc2('delz_inc',increments%delz_inc) + call interp_inc2('u_inc',increments%ua_inc) + call interp_inc2('v_inc',increments%va_inc) + do l=1,ntracers + call interp_inc2(trim(tracer_names(l))//'_inc',increments%tracer_inc(:,:,:,l)) + enddo + else + call interp_inc('T_inc',increments%temp_inc(:,:,:),jbeg,jend) + call interp_inc('delp_inc',increments%delp_inc(:,:,:),jbeg,jend) + call interp_inc('delz_inc',increments%delz_inc(:,:,:),jbeg,jend) + call interp_inc('u_inc',increments%ua_inc(:,:,:),jbeg,jend) ! can these be treated as scalars? + call interp_inc('v_inc',increments%va_inc(:,:,:),jbeg,jend) + do l=1,ntracers + call interp_inc(trim(tracer_names(l))//'_inc',increments%tracer_inc(:,:,:,l),jbeg,jend) + enddo + end if call close_ncfile(ncid) deallocate (wk3) @@ -517,6 +544,22 @@ subroutine interp_inc(field_name,var,jbeg,jend) enddo end subroutine interp_inc +subroutine interp_inc2(field_name,var) +! interpolate increment from GSI gaussian grid to cubed sphere +! everying is on the A-grid, earth relative + character(len=*), intent(in) :: field_name + real, dimension(is:ie,js:je,1:km), intent(inout) :: var + integer:: i1, i2, j1, k,j,i,ierr + call check_var_exists(ncid, field_name, ierr) + if (ierr == 0) then + call get_var3_r4( ncid, field_name, 1,im, js,je, 1,km, wk3 ) + else + if (is_master()) print *,'warning: no increment for ',trim(field_name),' found, assuming zero' + wk3 = 0. + end if + var(is:ie,js:je,1:km)=wk3(is:ie,:,:) +end subroutine interp_inc2 + end module fv_iau_mod From cf1544280704a3ac3a036a89c263fd08aeb59632 Mon Sep 17 00:00:00 2001 From: Dusan Jovic Date: Fri, 12 Jul 2024 13:37:27 +0000 Subject: [PATCH 18/29] Fix do loop bounds in a2b_ord2 routine in to avoid reading from uninitialized elements of qin array Input array qin in a2b_ord2 routine is declared as qin(is-ng:ie+ng,js-ng:je+ng) but only elements in (is-1:ie+1, js-1:je+1) range are initialized in a calling routine, for example pin array in adv_pe routine in dyn_core.F90 --- model/a2b_edge.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/model/a2b_edge.F90 b/model/a2b_edge.F90 index c4530a131..0c5de7eff 100644 --- a/model/a2b_edge.F90 +++ b/model/a2b_edge.F90 @@ -377,8 +377,8 @@ subroutine a2b_ord2(qin, qout, gridstruct, npx, npy, is, ie, js, je, ng, replace if (gridstruct%bounded_domain) then - do j=js-2,je+1+2 - do i=is-2,ie+1+2 + do j=js,je+1 + do i=is,ie+1 qout(i,j) = 0.25*(qin(i-1,j-1)+qin(i,j-1)+qin(i-1,j)+qin(i,j)) enddo enddo From 1a8f38b69953912385f502821c2e55345c3b3fd4 Mon Sep 17 00:00:00 2001 From: Dusan Jovic Date: Fri, 12 Jul 2024 13:48:15 +0000 Subject: [PATCH 19/29] Add '-init=snan,arrays' compiler flag to Intel debug Fortran flags (CMAKE_Fortran_FLAGS_DEBUG) --- cmake/compiler_flags_Intel_Fortran.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/compiler_flags_Intel_Fortran.cmake b/cmake/compiler_flags_Intel_Fortran.cmake index 0cfc70397..fbd667247 100644 --- a/cmake/compiler_flags_Intel_Fortran.cmake +++ b/cmake/compiler_flags_Intel_Fortran.cmake @@ -10,7 +10,7 @@ set(CMAKE_Fortran_FLAGS_REPRO "-O2 -debug minimal -fp-model consistent -qoverrid set(CMAKE_Fortran_FLAGS_RELEASE "-O2 -debug minimal -fp-model source -qoverride-limits -qopt-prefetch=3") -set(CMAKE_Fortran_FLAGS_DEBUG "-O0 -check -check noarg_temp_created -check nopointer -warn -warn noerrors -fp-stack-check -fstack-protector-all -fpe0 -debug -ftrapuv") +set(CMAKE_Fortran_FLAGS_DEBUG "-O0 -check -check noarg_temp_created -check nopointer -warn -warn noerrors -fp-stack-check -fstack-protector-all -fpe0 -debug -ftrapuv -init=snan,arrays") set(CMAKE_Fortran_LINK_FLAGS "") From 7bdc2c6f579e60f6417d525d58b47b2ec4b5f9bb Mon Sep 17 00:00:00 2001 From: "Bin.Liu" Date: Fri, 9 May 2025 09:20:56 -0500 Subject: [PATCH 20/29] 3D-TKE EMDF GFS PBL scheme related changes from FIU (Ping Zhu, Ping.Zhu@noaa.gov, and Kwun Y. Fung). --- driver/UFS/atmosphere.F90 | 65 +++++++- model/dyn_core.F90 | 328 +++++++++++++++++++++++++++++++++++++- model/fv_arrays.F90 | 25 ++- model/fv_control.F90 | 9 ++ model/fv_dynamics.F90 | 21 ++- model/nh_utils.F90 | 70 ++++++++ model/sw_core.F90 | 12 ++ 7 files changed, 523 insertions(+), 7 deletions(-) diff --git a/driver/UFS/atmosphere.F90 b/driver/UFS/atmosphere.F90 index b9abdb28d..eb8653149 100644 --- a/driver/UFS/atmosphere.F90 +++ b/driver/UFS/atmosphere.F90 @@ -182,6 +182,8 @@ module atmosphere_mod #ifdef GFS_TYPES use GFS_typedefs, only: IPD_control_type => GFS_control_type, kind_phys use GFS_typedefs, only: GFS_statein_type, GFS_stateout_type, GFS_sfcprop_type +! SA-3D-TKE (Samuel) +use GFS_typedefs, only: GFS_tbd_type #else use IPD_typedefs, only: IPD_data_type, IPD_control_type, kind_phys => IPD_kind_phys #endif @@ -693,11 +695,17 @@ subroutine atmosphere_dynamics ( Time ) ! Atm(n)%flagstruct%n_split, Atm(n)%flagstruct%q_split, & Atm(n)%u, Atm(n)%v, Atm(n)%w, Atm(n)%delz, & Atm(n)%flagstruct%hydrostatic, & +!The following variable is used for SA-3D-TKE + Atm(n)%flagstruct%sa3dtke_dyco, & Atm(n)%pt , Atm(n)%delp, Atm(n)%q, Atm(n)%ps, & Atm(n)%pe, Atm(n)%pk, Atm(n)%peln, & Atm(n)%pkz, Atm(n)%phis, Atm(n)%q_con, & Atm(n)%omga, Atm(n)%ua, Atm(n)%va, Atm(n)%uc, & - Atm(n)%vc, Atm(n)%ak, Atm(n)%bk, Atm(n)%mfx, & + Atm(n)%vc, & +!The following three variables are used for SA-3D-TKE + Atm(n)%deform_1,Atm(n)%deform_2, & + Atm(n)%deform_3,Atm(n)%dku3d_h,Atm(n)%dku3d_e, & + Atm(n)%ak, Atm(n)%bk, Atm(n)%mfx, & Atm(n)%mfy , Atm(n)%cx, Atm(n)%cy, Atm(n)%ze0, & Atm(n)%flagstruct%hybrid_z, & Atm(n)%gridstruct, Atm(n)%flagstruct, & @@ -1894,9 +1902,14 @@ subroutine adiabatic_init(zvir,nudge_dz,time) Atm(mygrid)%ptop, Atm(mygrid)%ks, nq, n_split_loc, & Atm(mygrid)%flagstruct%q_split, Atm(mygrid)%u, Atm(mygrid)%v, Atm(mygrid)%w, & Atm(mygrid)%delz, Atm(mygrid)%flagstruct%hydrostatic, & +!The following variable is used for SA-3D-TKE + Atm(mygrid)%flagstruct%sa3dtke_dyco, & Atm(mygrid)%pt, Atm(mygrid)%delp, Atm(mygrid)%q, Atm(mygrid)%ps, & Atm(mygrid)%pe, Atm(mygrid)%pk, Atm(mygrid)%peln, Atm(mygrid)%pkz, Atm(mygrid)%phis, & Atm(mygrid)%q_con, Atm(mygrid)%omga, Atm(mygrid)%ua, Atm(mygrid)%va, Atm(mygrid)%uc, Atm(mygrid)%vc, & +!The following three variables are used for SA-3D-TKE + Atm(mygrid)%deform_1,Atm(mygrid)%deform_2, & + Atm(mygrid)%deform_3,Atm(mygrid)%dku3d_h,Atm(mygrid)%dku3d_e, & Atm(mygrid)%ak, Atm(mygrid)%bk, Atm(mygrid)%mfx, Atm(mygrid)%mfy, & Atm(mygrid)%cx, Atm(mygrid)%cy, Atm(mygrid)%ze0, Atm(mygrid)%flagstruct%hybrid_z, & Atm(mygrid)%gridstruct, Atm(mygrid)%flagstruct, & @@ -1909,9 +1922,14 @@ subroutine adiabatic_init(zvir,nudge_dz,time) Atm(mygrid)%ptop, Atm(mygrid)%ks, nq, n_split_loc, & Atm(mygrid)%flagstruct%q_split, Atm(mygrid)%u, Atm(mygrid)%v, Atm(mygrid)%w, & Atm(mygrid)%delz, Atm(mygrid)%flagstruct%hydrostatic, & +!The following variable is used for SA-3D-TKE + Atm(mygrid)%flagstruct%sa3dtke_dyco, & Atm(mygrid)%pt, Atm(mygrid)%delp, Atm(mygrid)%q, Atm(mygrid)%ps, & Atm(mygrid)%pe, Atm(mygrid)%pk, Atm(mygrid)%peln, Atm(mygrid)%pkz, Atm(mygrid)%phis, & Atm(mygrid)%q_con, Atm(mygrid)%omga, Atm(mygrid)%ua, Atm(mygrid)%va, Atm(mygrid)%uc, Atm(mygrid)%vc, & +!The following three variables are used for SA-3D-TKE + Atm(mygrid)%deform_1, Atm(mygrid)%deform_2, & + Atm(mygrid)%deform_3,Atm(mygrid)%dku3d_h,Atm(mygrid)%dku3d_e, & Atm(mygrid)%ak, Atm(mygrid)%bk, Atm(mygrid)%mfx, Atm(mygrid)%mfy, & Atm(mygrid)%cx, Atm(mygrid)%cy, Atm(mygrid)%ze0, Atm(mygrid)%flagstruct%hybrid_z, & Atm(mygrid)%gridstruct, Atm(mygrid)%flagstruct, & @@ -1985,9 +2003,14 @@ subroutine adiabatic_init(zvir,nudge_dz,time) Atm(mygrid)%ptop, Atm(mygrid)%ks, nq, Atm(mygrid)%flagstruct%n_split, & Atm(mygrid)%flagstruct%q_split, Atm(mygrid)%u, Atm(mygrid)%v, Atm(mygrid)%w, & Atm(mygrid)%delz, Atm(mygrid)%flagstruct%hydrostatic, & +!The following variable is used for SA-3D-TKE + Atm(mygrid)%flagstruct%sa3dtke_dyco, & Atm(mygrid)%pt, Atm(mygrid)%delp, Atm(mygrid)%q, Atm(mygrid)%ps, & Atm(mygrid)%pe, Atm(mygrid)%pk, Atm(mygrid)%peln, Atm(mygrid)%pkz, Atm(mygrid)%phis, & Atm(mygrid)%q_con, Atm(mygrid)%omga, Atm(mygrid)%ua, Atm(mygrid)%va, Atm(mygrid)%uc, Atm(mygrid)%vc, & +!The following three variables are used for SA-3D-TKE + Atm(mygrid)%deform_1,Atm(mygrid)%deform_2, & + Atm(mygrid)%deform_3,Atm(mygrid)%dku3d_h,Atm(mygrid)%dku3d_e, & Atm(mygrid)%ak, Atm(mygrid)%bk, Atm(mygrid)%mfx, Atm(mygrid)%mfy, & Atm(mygrid)%cx, Atm(mygrid)%cy, Atm(mygrid)%ze0, Atm(mygrid)%flagstruct%hybrid_z, & Atm(mygrid)%gridstruct, Atm(mygrid)%flagstruct, & @@ -1999,9 +2022,14 @@ subroutine adiabatic_init(zvir,nudge_dz,time) Atm(mygrid)%ptop, Atm(mygrid)%ks, nq, Atm(mygrid)%flagstruct%n_split, & Atm(mygrid)%flagstruct%q_split, Atm(mygrid)%u, Atm(mygrid)%v, Atm(mygrid)%w, & Atm(mygrid)%delz, Atm(mygrid)%flagstruct%hydrostatic, & +!The following variable is used for SA-3D-TKE + Atm(mygrid)%flagstruct%sa3dtke_dyco, & Atm(mygrid)%pt, Atm(mygrid)%delp, Atm(mygrid)%q, Atm(mygrid)%ps, & Atm(mygrid)%pe, Atm(mygrid)%pk, Atm(mygrid)%peln, Atm(mygrid)%pkz, Atm(mygrid)%phis, & Atm(mygrid)%q_con, Atm(mygrid)%omga, Atm(mygrid)%ua, Atm(mygrid)%va, Atm(mygrid)%uc, Atm(mygrid)%vc, & +!The following three variables are used for SA-3D-TKE + Atm(mygrid)%deform_1,Atm(mygrid)%deform_2, & + Atm(mygrid)%deform_3,Atm(mygrid)%dku3d_h,Atm(mygrid)%dku3d_e, & Atm(mygrid)%ak, Atm(mygrid)%bk, Atm(mygrid)%mfx, Atm(mygrid)%mfy, & Atm(mygrid)%cx, Atm(mygrid)%cy, Atm(mygrid)%ze0, Atm(mygrid)%flagstruct%hybrid_z, & Atm(mygrid)%gridstruct, Atm(mygrid)%flagstruct, & @@ -2064,6 +2092,7 @@ end subroutine adiabatic_init !>@detail Performs a mass adjustment to be consistent with the !! GFS physics and if necessary, converts quantities to hydrostatic !! representation. +!! SA-3D-TKE (added IPD_Tbd in the input) (kyf) #if defined(OVERLOAD_R4) #define _DBL_(X) DBLE(X) #define _RL_(X) REAL(X,KIND=4) @@ -2071,9 +2100,11 @@ end subroutine adiabatic_init #define _DBL_(X) X #define _RL_(X) X #endif - subroutine atmos_phys_driver_statein (IPD_Control, IPD_Statein, Atm_block,flip_vc) + subroutine atmos_phys_driver_statein (IPD_Control, IPD_Statein, IPD_Tbd, Atm_block,flip_vc) type (IPD_control_type), intent(in) :: IPD_Control type (GFS_statein_type), intent(inout) :: IPD_Statein + ! SA-3D-TKE (added IPD_Tbd) (kyf) + type (GFS_tbd_type), intent(inout) :: IPD_Tbd type (block_control_type), intent(in) :: Atm_block logical, intent(in) :: flip_vc !-------------------------------------- @@ -2113,8 +2144,10 @@ subroutine atmos_phys_driver_statein (IPD_Control, IPD_Statein, Atm_block,flip_v !--------------------------------------------------------------------- ! use most up to date atmospheric properties when running serially !--------------------------------------------------------------------- +! SA-3D-TKE added IPD_Tbd (kyf) !$OMP parallel do default (none) & !$OMP shared (Atm_block, Atm, IPD_Control, IPD_Statein, npz, nq, ncnst, sphum, liq_wat, & +!$OMP IPD_Tbd, & !$OMP ice_wat, rainwat, snowwat, graupel, pk0inv, ptop, & !$OMP pktop, zvir, mygrid, dnats, nq_adv, flip_vc) & #ifdef MULTI_GASES @@ -2130,6 +2163,27 @@ subroutine atmos_phys_driver_statein (IPD_Control, IPD_Statein, Atm_block,flip_v ! log(pe) <-- prsik blen = Atm_block%blksz(nb) +!The following is for SA-3D-TKE + if(Atm(mygrid)%flagstruct%sa3dtke_dyco) then +!The following is for SA-3D-TKE (pass dku to dyn_core) + do k = 1, npz + kz = npz+1-k + if(flip_vc) then + k1 = kz ! flipping the index + else + k1 = k + endif + do ix = 1, blen + i = Atm_block%index(nb)%ii(ix) + j = Atm_block%index(nb)%jj(ix) + ! SA-3D-TKE (added im) (kyf) + im = IPD_control%chunk_begin(nb)+ix-1 + ! SA-3D-TKE (modified IPD_Tbd and ix into im) (kyf) + Atm(mygrid)%dku3d_h(i,j,k) = IPD_Tbd%dku3d_h(im,k1) + Atm(mygrid)%dku3d_e(i,j,k) = IPD_Tbd%dku3d_e(im,k1) + enddo + enddo + endif do k = 1, npz !Indices for FV's vertical coordinate, for which 1 = top @@ -2154,6 +2208,13 @@ subroutine atmos_phys_driver_statein (IPD_Control, IPD_Statein, Atm_block,flip_v endif IPD_Statein%vvl(im,k) = _DBL_(_RL_(Atm(mygrid)%omga(i,j,k1))) IPD_Statein%prsl(im,k) = _DBL_(_RL_(Atm(mygrid)%delp(i,j,k1))) ! Total mass +!The following is for SA-3D-TKE + if(Atm(mygrid)%flagstruct%sa3dtke_dyco) then + IPD_Statein%def_1(ix,k) = _DBL_(_RL_(Atm(mygrid)%deform_1(i,j,k1))) + IPD_Statein%def_2(ix,k) = _DBL_(_RL_(Atm(mygrid)%deform_2(i,j,k1))) + IPD_Statein%def_3(ix,k) = _DBL_(_RL_(Atm(mygrid)%deform_3(i,j,k1))) + endif + if (Atm(mygrid)%flagstruct%do_skeb) IPD_Statein%diss_est(im,k) = _DBL_(_RL_(Atm(mygrid)%diss_est(i,j,k1))) if(flip_vc) then diff --git a/model/dyn_core.F90 b/model/dyn_core.F90 index ce53b8467..a965826f4 100644 --- a/model/dyn_core.F90 +++ b/model/dyn_core.F90 @@ -116,6 +116,7 @@ module dyn_core_mod use sw_core_mod, only: c_sw, d_sw, d_md use a2b_edge_mod, only: a2b_ord2, a2b_ord4 use nh_core_mod, only: Riem_Solver3, Riem_Solver_C, update_dz_c, update_dz_d, nh_bc + use nh_utils_mod, only: edge_profile1 ! KGao: for dudz,dvdz,dwdz calculations use tp_core_mod, only: copy_corners use fv_timing_mod, only: timing_on, timing_off use fv_diagnostics_mod, only: prt_maxmin, fv_time, prt_mxm @@ -136,6 +137,9 @@ module dyn_core_mod use fv_regional_mod, only: regional_boundary_update use fv_regional_mod, only: current_time_in_seconds, bc_time_interval use fv_regional_mod, only: delz_regBC ! TEMPORARY --- lmh +!The following is for SA-3D-TKE + use tracer_manager_mod, only: get_tracer_names, get_number_tracers, get_tracer_index + use field_manager_mod, only: MODEL_ATMOS #ifdef SW_DYNAMICS use test_cases_mod, only: test_case, case9_forcing1, case9_forcing2 @@ -178,7 +182,11 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, #endif grav, hydrostatic, & u, v, w, delz, pt, q, delp, pe, pk, phis, ws, omga, ptop, pfull, ua, va, & - uc, vc, mfx, mfy, cx, cy, pkz, peln, q_con, ak, bk, & + uc, vc, & +!The following 6 variables are for SA-3D-TKE + sa3dtke_dyco, deform_1, deform_2, & + deform_3, dku3d_h, dku3d_e, & + mfx, mfy, cx, cy, pkz, peln, q_con, ak, bk, & ks, gridstruct, flagstruct, neststruct, idiag, bd, domain, & init_step, i_pack, end_step, diss_est,time_total) @@ -235,6 +243,14 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, real, intent(inout):: uc(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz) !< (uc, vc) are mostly used as the C grid winds real, intent(inout):: vc(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz) real, intent(inout), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: ua, va +!The following 6 variables are for SA-3D-TKE + logical, intent(IN) :: sa3dtke_dyco + real, intent(out), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: deform_1 + real, intent(out), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: deform_2 + real, intent(out), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: deform_3 + real, intent(inout), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: dku3d_h + real, intent(inout), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: dku3d_e + real, intent(inout):: q_con(bd%isd:, bd%jsd:, 1:) ! The Flux capacitors: accumulated Mass flux arrays @@ -288,6 +304,10 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, logical used real :: split_timestep_bc +!The following is for SA-3D-TKE + real:: tke(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + integer :: ntke + integer :: is, ie, js, je integer :: isd, ied, jsd, jed @@ -765,13 +785,16 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, endif +!The following is modified for SA-3D-TKE +!-- add sa3dtke_dyco and dku3d_h in parallel calculation call timing_on('d_sw') !$OMP parallel do default(none) shared(npz,flagstruct,nord_v,pfull,damp_vt,hydrostatic,last_step, & !$OMP is,ie,js,je,isd,ied,jsd,jed,omga,delp,gridstruct,npx,npy, & !$OMP ng,zh,vt,ptc,pt,u,v,w,uc,vc,ua,va,divgd,mfx,mfy,cx,cy, & +!$OMP sa3dtke_dyco, dku3d_h, & !$OMP crx,cry,xfx,yfx,q_con,zvir,sphum,nq,q,dt,bd,rdt,iep1,jep1, & -!$OMP heat_source,diss_est,ptop,first_call) & +!$OMP heat_source,diss_est,ptop,first_call) & !$OMP private(nord_k, nord_w, nord_t, damp_w, damp_t, d2_divg, & !$OMP d_con_k,kgb, hord_m, hord_v, hord_t, hord_p, wk, heat_s,diss_e, z_rat) do k=1,npz @@ -883,6 +906,8 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, u(isd,jsd,k), v(isd,jsd,k), w(isd:,jsd:,k), uc(isd,jsd,k), & vc(isd,jsd,k), ua(isd,jsd,k), va(isd,jsd,k), divgd(isd,jsd,k), & mfx(is, js, k), mfy(is, js, k), cx(is, jsd,k), cy(isd,js, k), & +!The following 2 variables are for SA-3D-TKE + sa3dtke_dyco, dku3d_h(isd, jsd, k), & crx(is, jsd,k), cry(isd,js, k), xfx(is, jsd,k), yfx(isd,js, k), & #ifdef USE_COND q_con(isd:,jsd:,k), z_rat(isd,jsd), & @@ -1370,6 +1395,24 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, enddo ! time split loop first_call=.false. !----------------------------------------------------- +!The following is for SA-3D-TKE +!--calculating shear deformation and TKE transport for 3d TKE scheme + if(sa3dtke_dyco) then + if ( end_step ) then + ntke = get_tracer_index(MODEL_ATMOS, 'sgs_tke') + call mpp_update_domains(q(:,:,:,ntke), domain) + call mpp_update_domains(dku3d_e(:,:,:), domain) ! update dku3d_e at halo + + call diff3d(npx, npy, npz, nq, ua, va, w, & + q, deform_1, deform_2, deform_3, dku3d_e, & + zh, dp_ref, gridstruct, bd) + call mpp_update_domains(deform_1, domain) + call mpp_update_domains(deform_2, domain) + call mpp_update_domains(deform_3, domain) + endif + endif + + if ( nq > 0 .and. .not. flagstruct%inline_q ) then call timing_on('COMM_TOTAL') call timing_on('COMM_TRACER') @@ -2112,6 +2155,287 @@ subroutine one_grad_p(u, v, pk, gz, divg2, delp, dt, ng, gridstruct, bd, npx, np end subroutine one_grad_p +!The following is for SA-3D-TKE +subroutine diff3d(npx, npy, npz, nq, ua, va, w, q, & + deform_1, deform_2, deform_3, dku3d_e, & + zh, dp_ref, gridstruct, bd) + + use fv_arrays_mod, only: fv_grid_type + + integer, intent(in) :: nq, npx, npy, npz + type(fv_grid_bounds_type), intent(IN) :: bd + real, intent(in) :: ua(bd%isd:bd%ied, bd%jsd:bd%jed, npz) + real, intent(in) :: va(bd%isd:bd%ied, bd%jsd:bd%jed, npz) + real, intent(in) :: w(bd%isd:bd%ied, bd%jsd:bd%jed, npz) + real, intent(in) :: q(bd%isd:bd%ied, bd%jsd:bd%jed, npz, nq) + real, intent(in) :: zh(bd%isd:bd%ied, bd%jsd:bd%jed, npz+1) + real, intent(in) :: dku3d_e(bd%isd:bd%ied ,bd%jsd:bd%jed,npz) + + !real, intent(in) :: gz(bd%is:,bd%js:,1:) ! KGao: dims may not be right; zh is now used + real, intent(in) :: dp_ref(npz) + + real, intent(out) :: deform_1(bd%isd:bd%ied, bd%jsd:bd%jed, npz) + real, intent(out) :: deform_2(bd%isd:bd%ied, bd%jsd:bd%jed, npz) + real, intent(out) :: deform_3(bd%isd:bd%ied, bd%jsd:bd%jed, npz) + + type(fv_grid_type), intent(IN), target :: gridstruct + +! local + + real, pointer, dimension(:,:) :: dx, dy, rarea + integer :: is, ie, js, je, isd, ied, jsd, jed, i, j, k + real:: dudx(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: dudy(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: dvdx(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: dvdy(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: dwdx(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: dwdy(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: dudz(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: dvdz(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: dwdz(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: ut(bd%isd:bd%ied+1,bd%jsd:bd%jed) + real:: vt(bd%isd:bd%ied, bd%jsd:bd%jed+1) + real:: u_e(bd%is:bd%ie,npz+1) + real:: v_e(bd%is:bd%ie,npz+1) + real:: w_e(bd%is:bd%ie,npz+1) + + real:: dkdx(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + real:: dkdy(bd%isd:bd%ied,bd%jsd:bd%jed,npz) + + integer :: ntke + real :: tke_1(bd%isd:bd%ied+2, bd%jsd:bd%jed+2) + real :: dedy_1(bd%isd:bd%ied+1, bd%jsd:bd%jed+1, npz) + real :: dedy_2(bd%isd:bd%ied, bd%jsd:bd%jed, npz) + real :: dedx_1(bd%isd:bd%ied+1, bd%jsd:bd%jed+1, npz) + real :: dedx_2(bd%isd:bd%ied, bd%jsd:bd%jed, npz) + real :: tke_2(npz) + real :: dedz_1(npz) + real :: dedz_2(bd%isd:bd%ied, bd%jsd:bd%jed,npz) + real :: tmp + + real :: dz + + dx => gridstruct%dx + dy => gridstruct%dy + rarea => gridstruct%rarea + + is = bd%is + ie = bd%ie + js = bd%js + je = bd%je + isd = bd%isd + ied = bd%ied + jsd = bd%jsd + jed = bd%jed + +!=========================================================== +! Calculate deform_1 and deform_2 +! deform_1 = 2*dw/dz**2+(du/dz**2+dv/dz**2) +! +(dw/dx*du/dz+dw/dy*dv/dz) +! deform_2 = 2*(du/dx**2+dv/dy**2)+(du/dy+dv/dx)**2 +! +(dw/dx**2+dw/dy**2)+(dw/dx*du/dz+dw/dy*dv/dz) +!=========================================================== + +! KGao: make ut and vt private as suggested by Lucas + +!$OMP parallel do default(none) shared(npz,is,ie,js,je,ua,va,w,dx,dy,rarea, & +!$OMP dudx,dudy,dvdx,dvdy,dwdx,dwdy) & +!$OMP private(ut,vt) + do k=1,npz +!------------------------------------- +! get du/dy and dv/dx + do j=js,je+1 + do i=is,ie + vt(i,j) = ua(i,j,k)*dx(i,j) + enddo + enddo + do j=js,je + do i=is,ie+1 + ut(i,j) = va(i,j,k)*dy(i,j) + enddo + enddo + do j=js,je + do i=is,ie + dudy(i,j,k) = rarea(i,j)*(vt(i,j+1)-vt(i,j)) + dvdx(i,j,k) = rarea(i,j)*(ut(i+1,j)-ut(i,j)) + enddo + enddo +!------------------------------------- +! get du/dx and dv/dy + do j=js,je + do i=is,ie+1 + ut(i,j) = ua(i,j,k)*dy(i,j) + enddo + enddo + do j=js,je+1 + do i=is,ie + vt(i,j) = va(i,j,k)*dx(i,j) + enddo + enddo + do j=js,je + do i=is,ie + dudx(i,j,k) = rarea(i,j)*(ut(i+1,j)-ut(i,j)) + dvdy(i,j,k) = rarea(i,j)*(vt(i,j+1)-vt(i,j)) + enddo + enddo +!------------------------------------- +! get dw/dx and dw/dy + do j=js,je + do i=is,ie+1 + ut(i,j) = w(i,j,k)*dy(i,j) + enddo + enddo + do j=js,je+1 + do i=is,ie + vt(i,j) = w(i,j,k)*dx(i,j) + enddo + enddo + do j=js,je + do i=is,ie + dwdx(i,j,k) = rarea(i,j)*(ut(i+1,j)-ut(i,j)) + dwdy(i,j,k) = rarea(i,j)*(vt(i,j+1)-vt(i,j)) + enddo + enddo + enddo !z loop + +!------------------------------------- +! get du/dz, dv/dz and dw/dz + +! KGao: use Lucas's method as in compute_dudz() +!$OMP parallel do default(none) shared(npz,is,ie,js,je,ua,va,w, & +!$OMP zh,dp_ref,dudz,dvdz,dwdz) & +!$OMP private(u_e,v_e,w_e,dz) + do j=js,je + call edge_profile1(ua(is:ie,j,:), u_e, is, ie, npz, dp_ref, 1) + call edge_profile1(va(is:ie,j,:), v_e, is, ie, npz, dp_ref, 1) + call edge_profile1(w(is:ie,j,:), w_e, is, ie, npz, dp_ref, 1) + do k=1,npz + do i=is,ie + dz = zh(i,j,k) - zh(i,j,k+1) + dudz(i,j,k) = (u_e(i,k)-u_e(i,k+1))/dz + dvdz(i,j,k) = (v_e(i,k)-v_e(i,k+1))/dz + dwdz(i,j,k) = (w_e(i,k)-w_e(i,k+1))/dz + enddo + enddo + enddo + +!------------------------------------- +! get deform_1 and deform_2 based on all terms + +!$OMP parallel do default(none) shared(npz,is,ie,js,je,deform_1,deform_2, & +!$OMP dudx,dudy,dudz,dvdx,dvdy,dvdz,dwdx,dwdy,dwdz) & +!$OMP private(tmp) + do k=1,npz + do j=js,je + do i=is,ie + tmp=dwdx(i,j,k)*dudz(i,j,k)+dwdy(i,j,k)*dvdz(i,j,k) + deform_1(i,j,k)= 2*dwdz(i,j,k)**2+ & + dudz(i,j,k)**2+dvdz(i,j,k)**2+tmp + deform_2(i,j,k)= & + 2*(dudx(i,j,k)**2+dvdy(i,j,k)**2)+ & + (dudy(i,j,k)+dvdx(i,j,k))**2+ & + dwdx(i,j,k)**2+dwdy(i,j,k)**2+tmp + enddo + enddo + enddo + +!=========================================================== +! Calculate deform_3 +! deform_3 = d(kqde/dx)/dx + d(kqde/dy)/dy +! where e is tke +!=========================================================== + + ntke = get_tracer_index(MODEL_ATMOS, 'sgs_tke') + +! KGao: make the 2d temporay arrays private - as suggested by Lucas + +!$OMP parallel do default(none) shared(npz,is,ie,js,je,ua,va,w,q,dx,dy,rarea, & +!$OMP dku3d_e,dkdx,dkdy,dedx_1,dedy_1,ntke,dedy_2,dedx_2) & +!$OMP private(ut,vt,tke_1) + do k=1,npz +!------------------------------------- + do j=js,je+2 + do i=is,ie+2 + tke_1(i,j)=q(i,j,k,ntke) + enddo + enddo + do j=js,je+2 + do i=is,ie + vt(i,j)=tke_1(i,j)*dx(i,j) + enddo + enddo + do j=js,je+1 + do i=is,ie + dedy_1(i,j,k)=rarea(i,j)*(vt(i,j+1)-vt(i,j)) + enddo + enddo + do j=js,je+1 + do i=is,ie + vt(i,j)=dedy_1(i,j,k)*dx(i,j) + enddo + enddo + do j=js,je + do i=is,ie + dedy_2(i,j,k)=rarea(i,j)*(vt(i,j+1)-vt(i,j)) + enddo + enddo +!------------------------------------- + do j=js,je + do i=is,ie+2 + ut(i,j)=tke_1(i,j)*dy(i,j) + enddo + enddo + do j=js,je + do i=is,ie+1 + dedx_1(i,j,k)=rarea(i,j)*(ut(i+1,j)-ut(i,j)) + enddo + enddo + do j=js,je + do i=is,ie+1 + ut(i,j)=dedx_1(i,j,k)*dy(i,j) + enddo + enddo + do j=js,je + do i=is,ie + dedx_2(i,j,k)=rarea(i,j)*(ut(i+1,j)-ut(i,j)) + enddo + enddo + do j=js,je+1 + do i=is,ie + vt(i,j)=dku3d_e(i,j,k)*dx(i,j) + enddo + enddo + do j=js,je + do i=is,ie + dkdy(i,j,k)=rarea(i,j)*(vt(i,j+1)-vt(i,j)) + enddo + enddo + do j=js,je + do i=is,ie+1 + ut(i,j)=dku3d_e(i,j,k)*dy(i,j) + enddo + enddo + do j=js,je + do i=is,ie + dkdx(i,j,k)=rarea(i,j)*(ut(i+1,j)-ut(i,j)) + enddo + enddo + enddo ! z loop +!------------------------------------- +! get deform_3 based on all terms + +!$OMP parallel do default(none) shared(npz,is,ie,js,je,deform_3, & +!$OMP dku3d_e,dkdx,dkdy,dedx_1,dedy_1,ntke,dedx_2,dedy_2) + do k=1,npz + do j=js,je + do i=is,ie + deform_3(i,j,k)=dku3d_e(i,j,k)*(dedx_2(i,j,k)+dedy_2(i,j,k)) & + +dkdx(i,j,k)*dedx_1(i,j,k)+dkdy(i,j,k)*dedy_1(i,j,k) + enddo + enddo + enddo + +end subroutine diff3d subroutine grad1_p_update(divg2, u, v, pk, gz, dt, ng, gridstruct, bd, npx, npy, npz, ptop, beta, a2b_ord) diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index 26a9cf153..969ddcf37 100644 --- a/model/fv_arrays.F90 +++ b/model/fv_arrays.F90 @@ -929,6 +929,8 @@ module fv_arrays_mod logical :: regional_bcs_from_gsi = .false. !< Default setting for writing restart files with boundary rows logical :: pass_full_omega_to_physics_in_non_hydrostatic_mode = .false. !< Default to passing local omega to physics in non-hydrostatic +!This logical variable is used for SA-3D-TKE + logical :: sa3dtke_dyco = .false. !>Convenience pointers integer, pointer :: grid_number @@ -1284,7 +1286,13 @@ module fv_arrays_mod real, _ALLOCATABLE :: va(:,:,:) _NULL real, _ALLOCATABLE :: uc(:,:,:) _NULL !< (uc, vc) are mostly used as the C grid winds real, _ALLOCATABLE :: vc(:,:,:) _NULL - +!3D-SA-TKE + real, _ALLOCATABLE :: deform_1(:,:,:) _NULL !< horizontal deformation + real, _ALLOCATABLE :: deform_2(:,:,:) _NULL !< vertical deformation + real, _ALLOCATABLE :: deform_3(:,:,:) _NULL !< 3D TKE transport & pressure correlation + real, _ALLOCATABLE :: dku3d_h(:,:,:) _NULL !< 3D Horizontal Eddy Diffusivity for Momentum + real, _ALLOCATABLE :: dku3d_e(:,:,:) _NULL !< 3D Eddy Diffusivity for TKE +!3D-SA-TKE-end real, _ALLOCATABLE :: ak(:) _NULL real, _ALLOCATABLE :: bk(:) _NULL @@ -1517,6 +1525,14 @@ subroutine allocate_fv_atmos_type(Atm, isd_in, ied_in, jsd_in, jed_in, is_in, ie allocate ( Atm%va(isd:ied ,jsd:jed ,npz) ) allocate ( Atm%uc(isd:ied+1,jsd:jed ,npz) ) allocate ( Atm%vc(isd:ied ,jsd:jed+1,npz) ) +!3D-SA-TKE + allocate ( Atm%deform_1(isd:ied ,jsd:jed ,npz) ) + allocate ( Atm%deform_2(isd:ied ,jsd:jed ,npz) ) + allocate ( Atm%deform_3(isd:ied ,jsd:jed ,npz) ) + allocate ( Atm%dku3d_h(isd:ied ,jsd:jed, npz) ) + allocate ( Atm%dku3d_e(isd:ied ,jsd:jed, npz) ) +!3D-SA-TKE-end + ! For tracer transport: allocate ( Atm%mfx(is:ie+1, js:je, npz) ) allocate ( Atm%mfy(is:ie , js:je+1,npz) ) @@ -1569,6 +1585,13 @@ subroutine allocate_fv_atmos_type(Atm, isd_in, ied_in, jsd_in, jed_in, is_in, ie Atm%va(i,j,k) = real_big Atm%pt(i,j,k) = real_big Atm%delp(i,j,k) = real_big +!3D-SA-TKE + Atm%deform_1(i,j,k) = 0. + Atm%deform_2(i,j,k) = 0. + Atm%deform_3(i,j,k) = 0. + Atm%dku3d_h(i,j,k) = 0. + Atm%dku3d_e(i,j,k) = 0. +!3D-SA-TKE-end #ifdef USE_COND Atm%q_con(i,j,k) = 0. #endif diff --git a/model/fv_control.F90 b/model/fv_control.F90 index 2ead0274c..f8eca5562 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -394,6 +394,9 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split, integer, pointer :: ndims +!The following is used for SA-3D-TKE + logical , pointer :: sa3dtke_dyco + real(kind=R_GRID), pointer :: dx_const real(kind=R_GRID), pointer :: dy_const real(kind=R_GRID), pointer :: deglon_start, deglon_stop, & ! boundaries of latlon patch @@ -952,6 +955,10 @@ subroutine set_namelist_pointers(Atm) external_eta => Atm%flagstruct%external_eta read_increment => Atm%flagstruct%read_increment increment_file_on_native_grid => Atm%flagstruct%increment_file_on_native_grid + +!The following is used for SA-3D-TKE + sa3dtke_dyco => Atm%flagstruct%sa3dtke_dyco + hydrostatic => Atm%flagstruct%hydrostatic phys_hydrostatic => Atm%flagstruct%phys_hydrostatic use_hydro_pressure => Atm%flagstruct%use_hydro_pressure @@ -1081,6 +1088,8 @@ subroutine read_namelist_fv_core_nml(Atm) consv_te, fill, filter_phys, fill_dp, fill_wz, fill_gfs, consv_am, RF_fast, & range_warn, dwind_2d, inline_q, z_tracer, reproduce_sum, adiabatic, do_vort_damp, no_dycore, & tau, tau_w, fast_tau_w_sec, tau_h2o, rf_cutoff, nf_omega, hydrostatic, fv_sg_adj, sg_cutoff, breed_vortex_inline, & +!The following is used for SA-3D-TKE + sa3dtke_dyco, & na_init, nudge_dz, hybrid_z, Make_NH, n_zs_filter, nord_zs_filter, full_zs_filter, reset_eta, & pnats, dnats, dnrts, a2b_ord, remap_t, p_ref, d2_bg_k1, d2_bg_k2, & c2l_ord, dx_const, dy_const, umax, deglat, & diff --git a/model/fv_dynamics.F90 b/model/fv_dynamics.F90 index 36b048136..2ad1de5b2 100644 --- a/model/fv_dynamics.F90 +++ b/model/fv_dynamics.F90 @@ -187,8 +187,13 @@ module fv_dynamics_mod subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, & reproduce_sum, kappa, cp_air, zvir, ptop, ks, ncnst, n_split, & - q_split, u, v, w, delz, hydrostatic, pt, delp, q, & + q_split, u, v, w, delz, hydrostatic, & +!The following variable is for SA-3D-TKE + sa3dtke_dyco, & + pt, delp, q, & ps, pe, pk, peln, pkz, phis, q_con, omga, ua, va, uc, vc, & +!The following three variables are for SA-3D-TKE + deform_1, deform_2, deform_3, dku3d_h, dku3d_e, & ak, bk, mfx, mfy, cx, cy, ze0, hybrid_z, & gridstruct, flagstruct, neststruct, idiag, bd, & parent_grid, domain, diss_est, inline_mp) @@ -255,6 +260,14 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, real, intent(inout) :: vc(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz) real, intent(inout), dimension(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz):: ua, va + !The following 4 variables are for SA-3D-TKE + logical, intent(in) :: sa3dtke_dyco + real, intent(out), dimension(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz):: deform_1 + real, intent(out), dimension(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz):: deform_2 + real, intent(out), dimension(bd%isd:bd%ied, bd%jsd:bd%jed, npz):: deform_3 + real, intent(inout), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: dku3d_h + real, intent(inout), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: dku3d_e + real, intent(in), dimension(npz+1):: ak, bk type(inline_mp_type), intent(inout) :: inline_mp @@ -666,7 +679,11 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, #endif grav, hydrostatic, & u, v, w, delz, pt, q, delp, pe, pk, phis, ws, omga, ptop, pfull, ua, va, & - uc, vc, mfx, mfy, cx, cy, pkz, peln, q_con, ak, bk, ks, & + uc, vc, & +!The following 6 variables are for SA-3D-TKE + sa3dtke_dyco, deform_1, deform_2, & + deform_3, dku3d_h, dku3d_e, & + mfx, mfy, cx, cy, pkz, peln, q_con, ak, bk, ks, & gridstruct, flagstruct, neststruct, idiag, bd, & domain, n_map==1, i_pack, GFDL_interstitial%last_step, diss_est,time_total) call timing_off('DYN_CORE') diff --git a/model/nh_utils.F90 b/model/nh_utils.F90 index 2b89aea48..7a42a2094 100644 --- a/model/nh_utils.F90 +++ b/model/nh_utils.F90 @@ -70,6 +70,9 @@ module nh_utils_mod public sim_solver, sim1_solver, sim3_solver public sim3p0_solver, rim_2d public Riem_Solver_c +!This is for sa3dtke + public edge_profile1 + real, parameter:: r3 = 1./3. @@ -1947,6 +1950,73 @@ subroutine edge_profile(q1, q2, q1e, q2e, i1, i2, j1, j2, j, km, dp0, uniform_gr end subroutine edge_profile + subroutine edge_profile1(q1, q1e, i1, i2, km, dp0, limiter) +! Edge profiles for a single scalar quantity + integer, intent(in):: i1, i2 + integer, intent(in):: km + integer, intent(in):: limiter + real, intent(in):: dp0(km) + real, intent(in), dimension(i1:i2,km):: q1 + real, intent(out), dimension(i1:i2,km+1):: q1e +!----------------------------------------------------------------------- + real, dimension(i1:i2,km+1):: qe1, gam ! edge values + real gak(km) + real bet, r2o3, r4o3 + real g0, gk, xt1, xt2, a_bot + integer i, k + +! Assuming grid varying in vertical only + g0 = dp0(2) / dp0(1) + xt1 = 2.*g0*(g0+1. ) + bet = g0*(g0+0.5) + do i=i1,i2 + qe1(i,1) = ( xt1*q1(i,1) + q1(i,2) ) / bet + gam(i,1) = ( 1. + g0*(g0+1.5) ) / bet + enddo + + do k=2,km + gk = dp0(k-1) / dp0(k) + do i=i1,i2 + bet = 2. + 2.*gk - gam(i,k-1) + qe1(i,k) = ( 3.*(q1(i,k-1)+gk*q1(i,k)) - qe1(i,k-1) ) / bet + gam(i,k) = gk / bet + enddo + enddo + + a_bot = 1. + gk*(gk+1.5) + xt1 = 2.*gk*(gk+1.) + do i=i1,i2 + xt2 = gk*(gk+0.5) - a_bot*gam(i,km) + qe1(i,km+1) = ( xt1*q1(i,km) + q1(i,km-1) - a_bot*qe1(i,km) ) / xt2 + enddo + + do k=km,1,-1 + do i=i1,i2 + qe1(i,k) = qe1(i,k) - gam(i,k)*qe1(i,k+1) + enddo + enddo + +!------------------ +! Apply constraints +!------------------ + if ( limiter/=0 ) then ! limit the top & bottom winds + do i=i1,i2 +! Top + if ( q1(i,1)*qe1(i,1) < 0. ) qe1(i,1) = 0. +! Surface: + if ( q1(i,km)*qe1(i,km+1) < 0. ) qe1(i,km+1) = 0. + enddo + endif + + do k=1,km+1 + do i=i1,i2 + q1e(i,k) = qe1(i,k) + enddo + enddo + + end subroutine edge_profile1 + + subroutine edge_profile_0grad(q1, q2, q1e, q2e, i1, i2, j1, j2, j, km, dp0, uniform_grid, limiter) ! Optimized for wind profile reconstruction: ! Added this option by Henry Juang and Xiaqiong Zhou 1/21/2021 diff --git a/model/sw_core.F90 b/model/sw_core.F90 index 6ed4a2e50..06291a07c 100644 --- a/model/sw_core.F90 +++ b/model/sw_core.F90 @@ -521,6 +521,8 @@ end subroutine c_sw !! and other prognostic varaiables. subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, & ua, va, divg_d, xflux, yflux, cx, cy, & +!The following 5 variables are for SA-3D-TKE + sa3dtke_dyco, dku3d_h, & crx_adv, cry_adv, xfx_adv, yfx_adv, q_con, z_rat, kgb, heat_source,diss_est, & zvir, sphum, nq, q, k, km, inline_q, & dt, hord_tr, hord_mt, hord_vt, hord_tm, hord_dp, nord, & @@ -537,6 +539,10 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, & real , intent(IN):: zvir real, intent(in):: damp_v, damp_w, damp_t, kgb type(fv_grid_bounds_type), intent(IN) :: bd +!The following 5 variables are for SA-3D-TKE + logical, intent(IN):: sa3dtke_dyco + real, intent(IN), dimension(bd%isd:bd%ied, bd%jsd:bd%jed):: dku3d_h + real, intent(inout):: divg_d(bd%isd:bd%ied+1,bd%jsd:bd%jed+1) !< divergence real, intent(IN), dimension(bd%isd:bd%ied, bd%jsd:bd%jed):: z_rat real, intent(INOUT), dimension(bd%isd:bd%ied, bd%jsd:bd%jed):: delp, pt, ua, va @@ -1453,7 +1459,13 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, & do j=js,je+1 do i=is,ie+1 +!The following is for SA-3D-TKE + if(sa3dtke_dyco) then + damp2 = abs(dt)*dku3d_h(i,j) + else damp2 = gridstruct%da_min_c*max(d2_bg, min(0.20, dddmp*vort(i,j))) ! del-2 + endif + vort(i,j) = damp2*delpc(i,j) + dd8*divg_d(i,j) ke(i,j) = ke(i,j) + vort(i,j) enddo From c50e3397fc5b74e30aca031e606396e4d6cbc9f0 Mon Sep 17 00:00:00 2001 From: "Kwun Y. Fung" Date: Fri, 23 May 2025 15:13:14 -0500 Subject: [PATCH 21/29] Update submodules fv3/atmos_cubed_sphere for 3dtke data structure --- driver/UFS/atmosphere.F90 | 41 ++++++++++++++---------------- model/dyn_core.F90 | 37 +++++++++++++--------------- model/fv_arrays.F90 | 52 +++++++++++++++++++++++---------------- model/fv_dynamics.F90 | 20 ++++++--------- model/sw_core.F90 | 10 +++++--- 5 files changed, 80 insertions(+), 80 deletions(-) diff --git a/driver/UFS/atmosphere.F90 b/driver/UFS/atmosphere.F90 index eb8653149..ac3524534 100644 --- a/driver/UFS/atmosphere.F90 +++ b/driver/UFS/atmosphere.F90 @@ -182,7 +182,7 @@ module atmosphere_mod #ifdef GFS_TYPES use GFS_typedefs, only: IPD_control_type => GFS_control_type, kind_phys use GFS_typedefs, only: GFS_statein_type, GFS_stateout_type, GFS_sfcprop_type -! SA-3D-TKE (Samuel) +! SA-3D-TKE (kyf) use GFS_typedefs, only: GFS_tbd_type #else use IPD_typedefs, only: IPD_data_type, IPD_control_type, kind_phys => IPD_kind_phys @@ -702,9 +702,8 @@ subroutine atmosphere_dynamics ( Time ) Atm(n)%pkz, Atm(n)%phis, Atm(n)%q_con, & Atm(n)%omga, Atm(n)%ua, Atm(n)%va, Atm(n)%uc, & Atm(n)%vc, & -!The following three variables are used for SA-3D-TKE - Atm(n)%deform_1,Atm(n)%deform_2, & - Atm(n)%deform_3,Atm(n)%dku3d_h,Atm(n)%dku3d_e, & +!The following variable is used for SA-3D-TKE (kyf) (modify for data structure) + Atm(n)%sa3dtke_var, & Atm(n)%ak, Atm(n)%bk, Atm(n)%mfx, & Atm(n)%mfy , Atm(n)%cx, Atm(n)%cy, Atm(n)%ze0, & Atm(n)%flagstruct%hybrid_z, & @@ -1907,9 +1906,8 @@ subroutine adiabatic_init(zvir,nudge_dz,time) Atm(mygrid)%pt, Atm(mygrid)%delp, Atm(mygrid)%q, Atm(mygrid)%ps, & Atm(mygrid)%pe, Atm(mygrid)%pk, Atm(mygrid)%peln, Atm(mygrid)%pkz, Atm(mygrid)%phis, & Atm(mygrid)%q_con, Atm(mygrid)%omga, Atm(mygrid)%ua, Atm(mygrid)%va, Atm(mygrid)%uc, Atm(mygrid)%vc, & -!The following three variables are used for SA-3D-TKE - Atm(mygrid)%deform_1,Atm(mygrid)%deform_2, & - Atm(mygrid)%deform_3,Atm(mygrid)%dku3d_h,Atm(mygrid)%dku3d_e, & +!The following variable is used for SA-3D-TKE (kyf) (modify for data structure) + Atm(mygrid)%sa3dtke_var, & Atm(mygrid)%ak, Atm(mygrid)%bk, Atm(mygrid)%mfx, Atm(mygrid)%mfy, & Atm(mygrid)%cx, Atm(mygrid)%cy, Atm(mygrid)%ze0, Atm(mygrid)%flagstruct%hybrid_z, & Atm(mygrid)%gridstruct, Atm(mygrid)%flagstruct, & @@ -1927,9 +1925,8 @@ subroutine adiabatic_init(zvir,nudge_dz,time) Atm(mygrid)%pt, Atm(mygrid)%delp, Atm(mygrid)%q, Atm(mygrid)%ps, & Atm(mygrid)%pe, Atm(mygrid)%pk, Atm(mygrid)%peln, Atm(mygrid)%pkz, Atm(mygrid)%phis, & Atm(mygrid)%q_con, Atm(mygrid)%omga, Atm(mygrid)%ua, Atm(mygrid)%va, Atm(mygrid)%uc, Atm(mygrid)%vc, & -!The following three variables are used for SA-3D-TKE - Atm(mygrid)%deform_1, Atm(mygrid)%deform_2, & - Atm(mygrid)%deform_3,Atm(mygrid)%dku3d_h,Atm(mygrid)%dku3d_e, & +!The following three variables are used for SA-3D-TKE (kyf) (modify for data structure) + Atm(mygrid)%sa3dtke_var, & Atm(mygrid)%ak, Atm(mygrid)%bk, Atm(mygrid)%mfx, Atm(mygrid)%mfy, & Atm(mygrid)%cx, Atm(mygrid)%cy, Atm(mygrid)%ze0, Atm(mygrid)%flagstruct%hybrid_z, & Atm(mygrid)%gridstruct, Atm(mygrid)%flagstruct, & @@ -2008,9 +2005,8 @@ subroutine adiabatic_init(zvir,nudge_dz,time) Atm(mygrid)%pt, Atm(mygrid)%delp, Atm(mygrid)%q, Atm(mygrid)%ps, & Atm(mygrid)%pe, Atm(mygrid)%pk, Atm(mygrid)%peln, Atm(mygrid)%pkz, Atm(mygrid)%phis, & Atm(mygrid)%q_con, Atm(mygrid)%omga, Atm(mygrid)%ua, Atm(mygrid)%va, Atm(mygrid)%uc, Atm(mygrid)%vc, & -!The following three variables are used for SA-3D-TKE - Atm(mygrid)%deform_1,Atm(mygrid)%deform_2, & - Atm(mygrid)%deform_3,Atm(mygrid)%dku3d_h,Atm(mygrid)%dku3d_e, & +!The following three variables are used for SA-3D-TKE (kyf) (modify for data structure) + Atm(mygrid)%sa3dtke_var, & Atm(mygrid)%ak, Atm(mygrid)%bk, Atm(mygrid)%mfx, Atm(mygrid)%mfy, & Atm(mygrid)%cx, Atm(mygrid)%cy, Atm(mygrid)%ze0, Atm(mygrid)%flagstruct%hybrid_z, & Atm(mygrid)%gridstruct, Atm(mygrid)%flagstruct, & @@ -2027,9 +2023,8 @@ subroutine adiabatic_init(zvir,nudge_dz,time) Atm(mygrid)%pt, Atm(mygrid)%delp, Atm(mygrid)%q, Atm(mygrid)%ps, & Atm(mygrid)%pe, Atm(mygrid)%pk, Atm(mygrid)%peln, Atm(mygrid)%pkz, Atm(mygrid)%phis, & Atm(mygrid)%q_con, Atm(mygrid)%omga, Atm(mygrid)%ua, Atm(mygrid)%va, Atm(mygrid)%uc, Atm(mygrid)%vc, & -!The following three variables are used for SA-3D-TKE - Atm(mygrid)%deform_1,Atm(mygrid)%deform_2, & - Atm(mygrid)%deform_3,Atm(mygrid)%dku3d_h,Atm(mygrid)%dku3d_e, & +!The following three variables are used for SA-3D-TKE (kyf) (modify for data structure) + Atm(mygrid)%sa3dtke_var, & Atm(mygrid)%ak, Atm(mygrid)%bk, Atm(mygrid)%mfx, Atm(mygrid)%mfy, & Atm(mygrid)%cx, Atm(mygrid)%cy, Atm(mygrid)%ze0, Atm(mygrid)%flagstruct%hybrid_z, & Atm(mygrid)%gridstruct, Atm(mygrid)%flagstruct, & @@ -2178,9 +2173,9 @@ subroutine atmos_phys_driver_statein (IPD_Control, IPD_Statein, IPD_Tbd, Atm_blo j = Atm_block%index(nb)%jj(ix) ! SA-3D-TKE (added im) (kyf) im = IPD_control%chunk_begin(nb)+ix-1 - ! SA-3D-TKE (modified IPD_Tbd and ix into im) (kyf) - Atm(mygrid)%dku3d_h(i,j,k) = IPD_Tbd%dku3d_h(im,k1) - Atm(mygrid)%dku3d_e(i,j,k) = IPD_Tbd%dku3d_e(im,k1) + ! SA-3D-TKE (modified IPD_Tbd and ix into im) (kyf) (modify for data structure) + Atm(mygrid)%sa3dtke_var%dku3d_h(i,j,k) = IPD_Tbd%dku3d_h(im,k1) + Atm(mygrid)%sa3dtke_var%dku3d_e(i,j,k) = IPD_Tbd%dku3d_e(im,k1) enddo enddo endif @@ -2208,11 +2203,11 @@ subroutine atmos_phys_driver_statein (IPD_Control, IPD_Statein, IPD_Tbd, Atm_blo endif IPD_Statein%vvl(im,k) = _DBL_(_RL_(Atm(mygrid)%omga(i,j,k1))) IPD_Statein%prsl(im,k) = _DBL_(_RL_(Atm(mygrid)%delp(i,j,k1))) ! Total mass -!The following is for SA-3D-TKE +!The following is for SA-3D-TKE (kyf) (modify for data structure) if(Atm(mygrid)%flagstruct%sa3dtke_dyco) then - IPD_Statein%def_1(ix,k) = _DBL_(_RL_(Atm(mygrid)%deform_1(i,j,k1))) - IPD_Statein%def_2(ix,k) = _DBL_(_RL_(Atm(mygrid)%deform_2(i,j,k1))) - IPD_Statein%def_3(ix,k) = _DBL_(_RL_(Atm(mygrid)%deform_3(i,j,k1))) + IPD_Statein%def_1(ix,k) = _DBL_(_RL_(Atm(mygrid)%sa3dtke_var%deform_1(i,j,k1))) + IPD_Statein%def_2(ix,k) = _DBL_(_RL_(Atm(mygrid)%sa3dtke_var%deform_2(i,j,k1))) + IPD_Statein%def_3(ix,k) = _DBL_(_RL_(Atm(mygrid)%sa3dtke_var%deform_3(i,j,k1))) endif if (Atm(mygrid)%flagstruct%do_skeb) IPD_Statein%diss_est(im,k) = _DBL_(_RL_(Atm(mygrid)%diss_est(i,j,k1))) diff --git a/model/dyn_core.F90 b/model/dyn_core.F90 index a965826f4..38d80de08 100644 --- a/model/dyn_core.F90 +++ b/model/dyn_core.F90 @@ -132,6 +132,7 @@ module dyn_core_mod use diag_manager_mod, only: send_data use fv_arrays_mod, only: fv_grid_type, fv_flags_type, fv_nest_type, fv_diag_type, & fv_grid_bounds_type, R_GRID, fv_nest_BC_type_3d + use fv_arrays_mod, only: sa3dtke_type ! for SA-3D-TKE (kyf) (modify for data structure) use boundary_mod, only: extrapolation_BC, nested_grid_BC_apply_intT use fv_regional_mod, only: regional_boundary_update @@ -183,9 +184,8 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, grav, hydrostatic, & u, v, w, delz, pt, q, delp, pe, pk, phis, ws, omga, ptop, pfull, ua, va, & uc, vc, & -!The following 6 variables are for SA-3D-TKE - sa3dtke_dyco, deform_1, deform_2, & - deform_3, dku3d_h, dku3d_e, & +!The following 2 variables are for SA-3D-TKE (kyf) (modify for data structure) + sa3dtke_dyco, sa3dtke_var, & mfx, mfy, cx, cy, pkz, peln, q_con, ak, bk, & ks, gridstruct, flagstruct, neststruct, idiag, bd, domain, & init_step, i_pack, end_step, diss_est,time_total) @@ -243,13 +243,9 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, real, intent(inout):: uc(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz) !< (uc, vc) are mostly used as the C grid winds real, intent(inout):: vc(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz) real, intent(inout), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: ua, va -!The following 6 variables are for SA-3D-TKE +!The following 2 variables are for SA-3D-TKE (kyf) (modify for data structure) logical, intent(IN) :: sa3dtke_dyco - real, intent(out), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: deform_1 - real, intent(out), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: deform_2 - real, intent(out), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: deform_3 - real, intent(inout), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: dku3d_h - real, intent(inout), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: dku3d_e + type(sa3dtke_type), intent(inout) :: sa3dtke_var real, intent(inout):: q_con(bd%isd:, bd%jsd:, 1:) @@ -787,12 +783,13 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, !The following is modified for SA-3D-TKE !-- add sa3dtke_dyco and dku3d_h in parallel calculation +! replce dku3d_h by sa3dtke_var (kyf) (modify for data structure) call timing_on('d_sw') !$OMP parallel do default(none) shared(npz,flagstruct,nord_v,pfull,damp_vt,hydrostatic,last_step, & !$OMP is,ie,js,je,isd,ied,jsd,jed,omga,delp,gridstruct,npx,npy, & !$OMP ng,zh,vt,ptc,pt,u,v,w,uc,vc,ua,va,divgd,mfx,mfy,cx,cy, & -!$OMP sa3dtke_dyco, dku3d_h, & +!$OMP sa3dtke_dyco, sa3dtke_var, & !$OMP crx,cry,xfx,yfx,q_con,zvir,sphum,nq,q,dt,bd,rdt,iep1,jep1, & !$OMP heat_source,diss_est,ptop,first_call) & !$OMP private(nord_k, nord_w, nord_t, damp_w, damp_t, d2_divg, & @@ -906,8 +903,8 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, u(isd,jsd,k), v(isd,jsd,k), w(isd:,jsd:,k), uc(isd,jsd,k), & vc(isd,jsd,k), ua(isd,jsd,k), va(isd,jsd,k), divgd(isd,jsd,k), & mfx(is, js, k), mfy(is, js, k), cx(is, jsd,k), cy(isd,js, k), & -!The following 2 variables are for SA-3D-TKE - sa3dtke_dyco, dku3d_h(isd, jsd, k), & +!The following 2 variables are for SA-3D-TKE (kyf) (modify for data structure) + sa3dtke_dyco, sa3dtke_var%dku3d_h(isd, jsd, k), & crx(is, jsd,k), cry(isd,js, k), xfx(is, jsd,k), yfx(isd,js, k), & #ifdef USE_COND q_con(isd:,jsd:,k), z_rat(isd,jsd), & @@ -1395,20 +1392,21 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, enddo ! time split loop first_call=.false. !----------------------------------------------------- -!The following is for SA-3D-TKE +!The following is for SA-3D-TKE (kyf) (modify for data structure) !--calculating shear deformation and TKE transport for 3d TKE scheme if(sa3dtke_dyco) then if ( end_step ) then ntke = get_tracer_index(MODEL_ATMOS, 'sgs_tke') call mpp_update_domains(q(:,:,:,ntke), domain) - call mpp_update_domains(dku3d_e(:,:,:), domain) ! update dku3d_e at halo + call mpp_update_domains(sa3dtke_var%dku3d_e(:,:,:), domain) ! update dku3d_e at halo call diff3d(npx, npy, npz, nq, ua, va, w, & - q, deform_1, deform_2, deform_3, dku3d_e, & + q, sa3dtke_var%deform_1, sa3dtke_var%deform_2, & + sa3dtke_var%deform_3, sa3dtke_var%dku3d_e, & zh, dp_ref, gridstruct, bd) - call mpp_update_domains(deform_1, domain) - call mpp_update_domains(deform_2, domain) - call mpp_update_domains(deform_3, domain) + call mpp_update_domains(sa3dtke_var%deform_1, domain) + call mpp_update_domains(sa3dtke_var%deform_2, domain) + call mpp_update_domains(sa3dtke_var%deform_3, domain) endif endif @@ -2177,7 +2175,6 @@ subroutine diff3d(npx, npy, npz, nq, ua, va, w, q, & real, intent(out) :: deform_1(bd%isd:bd%ied, bd%jsd:bd%jed, npz) real, intent(out) :: deform_2(bd%isd:bd%ied, bd%jsd:bd%jed, npz) real, intent(out) :: deform_3(bd%isd:bd%ied, bd%jsd:bd%jed, npz) - type(fv_grid_type), intent(IN), target :: gridstruct ! local @@ -2214,7 +2211,7 @@ subroutine diff3d(npx, npy, npz, nq, ua, va, w, q, & real :: tmp real :: dz - + dx => gridstruct%dx dy => gridstruct%dy rarea => gridstruct%rarea diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index 969ddcf37..a9946e8b9 100644 --- a/model/fv_arrays.F90 +++ b/model/fv_arrays.F90 @@ -1062,6 +1062,16 @@ module fv_arrays_mod logical :: BCfile_ne_is_open=.false. logical :: BCfile_sw_is_open=.false. end type fv_nest_type + + !3D-SA-TKE (kyf) (modify for data structure) + type sa3dtke_type + real, _ALLOCATABLE :: deform_1(:,:,:) _NULL !< horizontal deformation + real, _ALLOCATABLE :: deform_2(:,:,:) _NULL !< vertical deformation + real, _ALLOCATABLE :: deform_3(:,:,:) _NULL !< 3D TKE transport & pressure correlation + real, _ALLOCATABLE :: dku3d_h(:,:,:) _NULL !< 3D Horizontal Eddy Diffusivity for Momentum + real, _ALLOCATABLE :: dku3d_e(:,:,:) _NULL !< 3D Eddy Diffusivity for TKE + end type sa3dtke_type + !3D-SA-TKE-end type inline_mp_type real, _ALLOCATABLE :: prer(:,:) _NULL @@ -1286,13 +1296,6 @@ module fv_arrays_mod real, _ALLOCATABLE :: va(:,:,:) _NULL real, _ALLOCATABLE :: uc(:,:,:) _NULL !< (uc, vc) are mostly used as the C grid winds real, _ALLOCATABLE :: vc(:,:,:) _NULL -!3D-SA-TKE - real, _ALLOCATABLE :: deform_1(:,:,:) _NULL !< horizontal deformation - real, _ALLOCATABLE :: deform_2(:,:,:) _NULL !< vertical deformation - real, _ALLOCATABLE :: deform_3(:,:,:) _NULL !< 3D TKE transport & pressure correlation - real, _ALLOCATABLE :: dku3d_h(:,:,:) _NULL !< 3D Horizontal Eddy Diffusivity for Momentum - real, _ALLOCATABLE :: dku3d_e(:,:,:) _NULL !< 3D Eddy Diffusivity for TKE -!3D-SA-TKE-end real, _ALLOCATABLE :: ak(:) _NULL real, _ALLOCATABLE :: bk(:) _NULL @@ -1379,6 +1382,7 @@ module fv_arrays_mod integer :: atmos_axes(4) + type(sa3dtke_type) :: sa3dtke_var ! SA-3D TKE (kyf) (modify for data structure) type(inline_mp_type) :: inline_mp type(phys_diag_type) :: phys_diag type(nudge_diag_type) :: nudge_diag @@ -1525,13 +1529,16 @@ subroutine allocate_fv_atmos_type(Atm, isd_in, ied_in, jsd_in, jed_in, is_in, ie allocate ( Atm%va(isd:ied ,jsd:jed ,npz) ) allocate ( Atm%uc(isd:ied+1,jsd:jed ,npz) ) allocate ( Atm%vc(isd:ied ,jsd:jed+1,npz) ) -!3D-SA-TKE - allocate ( Atm%deform_1(isd:ied ,jsd:jed ,npz) ) - allocate ( Atm%deform_2(isd:ied ,jsd:jed ,npz) ) - allocate ( Atm%deform_3(isd:ied ,jsd:jed ,npz) ) - allocate ( Atm%dku3d_h(isd:ied ,jsd:jed, npz) ) - allocate ( Atm%dku3d_e(isd:ied ,jsd:jed, npz) ) -!3D-SA-TKE-end + + !3D-SA-TKE (kyf) (modify for data structure) + if ( Atm%flagstruct%sa3dtke_dyco ) then + allocate ( Atm%sa3dtke_var%deform_1(isd:ied ,jsd:jed ,npz) ) + allocate ( Atm%sa3dtke_var%deform_2(isd:ied ,jsd:jed ,npz) ) + allocate ( Atm%sa3dtke_var%deform_3(isd:ied ,jsd:jed ,npz) ) + allocate ( Atm%sa3dtke_var%dku3d_h(isd:ied ,jsd:jed, npz) ) + allocate ( Atm%sa3dtke_var%dku3d_e(isd:ied ,jsd:jed, npz) ) + endif + !3D-SA-TKE-end ! For tracer transport: allocate ( Atm%mfx(is:ie+1, js:je, npz) ) @@ -1585,13 +1592,16 @@ subroutine allocate_fv_atmos_type(Atm, isd_in, ied_in, jsd_in, jed_in, is_in, ie Atm%va(i,j,k) = real_big Atm%pt(i,j,k) = real_big Atm%delp(i,j,k) = real_big -!3D-SA-TKE - Atm%deform_1(i,j,k) = 0. - Atm%deform_2(i,j,k) = 0. - Atm%deform_3(i,j,k) = 0. - Atm%dku3d_h(i,j,k) = 0. - Atm%dku3d_e(i,j,k) = 0. -!3D-SA-TKE-end + + !3D-SA-TKE (kyf) (modify for data structure) + if ( Atm%flagstruct%sa3dtke_dyco ) then + Atm%sa3dtke_var%deform_1(i,j,k) = 0. + Atm%sa3dtke_var%deform_2(i,j,k) = 0. + Atm%sa3dtke_var%deform_3(i,j,k) = 0. + Atm%sa3dtke_var%dku3d_h(i,j,k) = 0. + Atm%sa3dtke_var%dku3d_e(i,j,k) = 0. + endif + !3D-SA-TKE-end #ifdef USE_COND Atm%q_con(i,j,k) = 0. #endif diff --git a/model/fv_dynamics.F90 b/model/fv_dynamics.F90 index 2ad1de5b2..fc98e9b6c 100644 --- a/model/fv_dynamics.F90 +++ b/model/fv_dynamics.F90 @@ -156,6 +156,7 @@ module fv_dynamics_mod use fv_regional_mod, only: current_time_in_seconds use boundary_mod, only: nested_grid_BC_apply_intT use fv_arrays_mod, only: fv_grid_type, fv_flags_type, fv_atmos_type, fv_nest_type, fv_diag_type, fv_grid_bounds_type, inline_mp_type + use fv_arrays_mod, only: sa3dtke_type ! for SA-3D-TKE (kyf) (modify for data structure) use fv_nwp_nudge_mod, only: do_adiabatic_init use time_manager_mod, only: get_time @@ -192,8 +193,8 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, sa3dtke_dyco, & pt, delp, q, & ps, pe, pk, peln, pkz, phis, q_con, omga, ua, va, uc, vc, & -!The following three variables are for SA-3D-TKE - deform_1, deform_2, deform_3, dku3d_h, dku3d_e, & +!The following variable is for SA-3D-TKE (kyf) (modify for data structure) + sa3dtke_var, & ak, bk, mfx, mfy, cx, cy, ze0, hybrid_z, & gridstruct, flagstruct, neststruct, idiag, bd, & parent_grid, domain, diss_est, inline_mp) @@ -260,16 +261,12 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, real, intent(inout) :: vc(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz) real, intent(inout), dimension(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz):: ua, va - !The following 4 variables are for SA-3D-TKE + !The following variable is for SA-3D-TKE (kyf) (modify for data structure) logical, intent(in) :: sa3dtke_dyco - real, intent(out), dimension(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz):: deform_1 - real, intent(out), dimension(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz):: deform_2 - real, intent(out), dimension(bd%isd:bd%ied, bd%jsd:bd%jed, npz):: deform_3 - real, intent(inout), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: dku3d_h - real, intent(inout), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: dku3d_e real, intent(in), dimension(npz+1):: ak, bk - + ! For SA-3D-TKE (kyf) (modify for data structure) + type(sa3dtke_type), intent(inout) :: sa3dtke_var type(inline_mp_type), intent(inout) :: inline_mp ! Accumulated Mass flux arrays: the "Flux Capacitor" @@ -680,9 +677,8 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, grav, hydrostatic, & u, v, w, delz, pt, q, delp, pe, pk, phis, ws, omga, ptop, pfull, ua, va, & uc, vc, & -!The following 6 variables are for SA-3D-TKE - sa3dtke_dyco, deform_1, deform_2, & - deform_3, dku3d_h, dku3d_e, & +!The following 2 variables are for SA-3D-TKE (kyf) (modify for data structure) + sa3dtke_dyco, sa3dtke_var, & mfx, mfy, cx, cy, pkz, peln, q_con, ak, bk, ks, & gridstruct, flagstruct, neststruct, idiag, bd, & domain, n_map==1, i_pack, GFDL_interstitial%last_step, diss_est,time_total) diff --git a/model/sw_core.F90 b/model/sw_core.F90 index 06291a07c..8e3bf0c85 100644 --- a/model/sw_core.F90 +++ b/model/sw_core.F90 @@ -50,6 +50,7 @@ module sw_core_mod deln_flux_explm, deln_flux_explm_udvd use fv_mp_mod, only: is_master, fill_corners, XDir, YDir use fv_arrays_mod, only: fv_grid_type, fv_grid_bounds_type, fv_flags_type + !use fv_arrays_mod, only: sa3dtke_type ! for SA-3D-TKE (kyf) (modify for data structure) use a2b_edge_mod, only: a2b_ord4 #ifdef SW_DYNAMICS @@ -521,7 +522,7 @@ end subroutine c_sw !! and other prognostic varaiables. subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, & ua, va, divg_d, xflux, yflux, cx, cy, & -!The following 5 variables are for SA-3D-TKE +!The following 2 variables are for SA-3D-TKE (kyf) (modify for data structure) sa3dtke_dyco, dku3d_h, & crx_adv, cry_adv, xfx_adv, yfx_adv, q_con, z_rat, kgb, heat_source,diss_est, & zvir, sphum, nq, q, k, km, inline_q, & @@ -539,9 +540,10 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, & real , intent(IN):: zvir real, intent(in):: damp_v, damp_w, damp_t, kgb type(fv_grid_bounds_type), intent(IN) :: bd -!The following 5 variables are for SA-3D-TKE +!The following 2 variables are for SA-3D-TKE (kyf) (modify for data structure) logical, intent(IN):: sa3dtke_dyco - real, intent(IN), dimension(bd%isd:bd%ied, bd%jsd:bd%jed):: dku3d_h + !type(sa3dtke_type), intent(in) :: sa3dtke_var + real, intent(IN), dimension(bd%isd:bd%ied, bd%jsd:bd%jed):: dku3d_h real, intent(inout):: divg_d(bd%isd:bd%ied+1,bd%jsd:bd%jed+1) !< divergence real, intent(IN), dimension(bd%isd:bd%ied, bd%jsd:bd%jed):: z_rat @@ -1461,7 +1463,7 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, & do i=is,ie+1 !The following is for SA-3D-TKE if(sa3dtke_dyco) then - damp2 = abs(dt)*dku3d_h(i,j) + damp2 = abs(dt)*dku3d_h(i,j) else damp2 = gridstruct%da_min_c*max(d2_bg, min(0.20, dddmp*vort(i,j))) ! del-2 endif From bd0e4141fccd40335c133ebe39768834b8056544 Mon Sep 17 00:00:00 2001 From: "Kwun Y. Fung" Date: Tue, 3 Jun 2025 00:02:14 -0500 Subject: [PATCH 22/29] Update based on stylistic comments and mpi use in dyn_core --- driver/UFS/atmosphere.F90 | 10 ---------- model/dyn_core.F90 | 42 +++++++++++++++++++-------------------- model/fv_dynamics.F90 | 7 ++----- model/sw_core.F90 | 5 ++--- 4 files changed, 24 insertions(+), 40 deletions(-) diff --git a/driver/UFS/atmosphere.F90 b/driver/UFS/atmosphere.F90 index ac3524534..61eb52d28 100644 --- a/driver/UFS/atmosphere.F90 +++ b/driver/UFS/atmosphere.F90 @@ -695,8 +695,6 @@ subroutine atmosphere_dynamics ( Time ) ! Atm(n)%flagstruct%n_split, Atm(n)%flagstruct%q_split, & Atm(n)%u, Atm(n)%v, Atm(n)%w, Atm(n)%delz, & Atm(n)%flagstruct%hydrostatic, & -!The following variable is used for SA-3D-TKE - Atm(n)%flagstruct%sa3dtke_dyco, & Atm(n)%pt , Atm(n)%delp, Atm(n)%q, Atm(n)%ps, & Atm(n)%pe, Atm(n)%pk, Atm(n)%peln, & Atm(n)%pkz, Atm(n)%phis, Atm(n)%q_con, & @@ -1901,8 +1899,6 @@ subroutine adiabatic_init(zvir,nudge_dz,time) Atm(mygrid)%ptop, Atm(mygrid)%ks, nq, n_split_loc, & Atm(mygrid)%flagstruct%q_split, Atm(mygrid)%u, Atm(mygrid)%v, Atm(mygrid)%w, & Atm(mygrid)%delz, Atm(mygrid)%flagstruct%hydrostatic, & -!The following variable is used for SA-3D-TKE - Atm(mygrid)%flagstruct%sa3dtke_dyco, & Atm(mygrid)%pt, Atm(mygrid)%delp, Atm(mygrid)%q, Atm(mygrid)%ps, & Atm(mygrid)%pe, Atm(mygrid)%pk, Atm(mygrid)%peln, Atm(mygrid)%pkz, Atm(mygrid)%phis, & Atm(mygrid)%q_con, Atm(mygrid)%omga, Atm(mygrid)%ua, Atm(mygrid)%va, Atm(mygrid)%uc, Atm(mygrid)%vc, & @@ -1920,8 +1916,6 @@ subroutine adiabatic_init(zvir,nudge_dz,time) Atm(mygrid)%ptop, Atm(mygrid)%ks, nq, n_split_loc, & Atm(mygrid)%flagstruct%q_split, Atm(mygrid)%u, Atm(mygrid)%v, Atm(mygrid)%w, & Atm(mygrid)%delz, Atm(mygrid)%flagstruct%hydrostatic, & -!The following variable is used for SA-3D-TKE - Atm(mygrid)%flagstruct%sa3dtke_dyco, & Atm(mygrid)%pt, Atm(mygrid)%delp, Atm(mygrid)%q, Atm(mygrid)%ps, & Atm(mygrid)%pe, Atm(mygrid)%pk, Atm(mygrid)%peln, Atm(mygrid)%pkz, Atm(mygrid)%phis, & Atm(mygrid)%q_con, Atm(mygrid)%omga, Atm(mygrid)%ua, Atm(mygrid)%va, Atm(mygrid)%uc, Atm(mygrid)%vc, & @@ -2000,8 +1994,6 @@ subroutine adiabatic_init(zvir,nudge_dz,time) Atm(mygrid)%ptop, Atm(mygrid)%ks, nq, Atm(mygrid)%flagstruct%n_split, & Atm(mygrid)%flagstruct%q_split, Atm(mygrid)%u, Atm(mygrid)%v, Atm(mygrid)%w, & Atm(mygrid)%delz, Atm(mygrid)%flagstruct%hydrostatic, & -!The following variable is used for SA-3D-TKE - Atm(mygrid)%flagstruct%sa3dtke_dyco, & Atm(mygrid)%pt, Atm(mygrid)%delp, Atm(mygrid)%q, Atm(mygrid)%ps, & Atm(mygrid)%pe, Atm(mygrid)%pk, Atm(mygrid)%peln, Atm(mygrid)%pkz, Atm(mygrid)%phis, & Atm(mygrid)%q_con, Atm(mygrid)%omga, Atm(mygrid)%ua, Atm(mygrid)%va, Atm(mygrid)%uc, Atm(mygrid)%vc, & @@ -2018,8 +2010,6 @@ subroutine adiabatic_init(zvir,nudge_dz,time) Atm(mygrid)%ptop, Atm(mygrid)%ks, nq, Atm(mygrid)%flagstruct%n_split, & Atm(mygrid)%flagstruct%q_split, Atm(mygrid)%u, Atm(mygrid)%v, Atm(mygrid)%w, & Atm(mygrid)%delz, Atm(mygrid)%flagstruct%hydrostatic, & -!The following variable is used for SA-3D-TKE - Atm(mygrid)%flagstruct%sa3dtke_dyco, & Atm(mygrid)%pt, Atm(mygrid)%delp, Atm(mygrid)%q, Atm(mygrid)%ps, & Atm(mygrid)%pe, Atm(mygrid)%pk, Atm(mygrid)%peln, Atm(mygrid)%pkz, Atm(mygrid)%phis, & Atm(mygrid)%q_con, Atm(mygrid)%omga, Atm(mygrid)%ua, Atm(mygrid)%va, Atm(mygrid)%uc, Atm(mygrid)%vc, & diff --git a/model/dyn_core.F90 b/model/dyn_core.F90 index 38d80de08..82ef021ca 100644 --- a/model/dyn_core.F90 +++ b/model/dyn_core.F90 @@ -184,8 +184,8 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, grav, hydrostatic, & u, v, w, delz, pt, q, delp, pe, pk, phis, ws, omga, ptop, pfull, ua, va, & uc, vc, & -!The following 2 variables are for SA-3D-TKE (kyf) (modify for data structure) - sa3dtke_dyco, sa3dtke_var, & +!The following variable is for SA-3D-TKE (kyf) (modify for data structure) + sa3dtke_var, & mfx, mfy, cx, cy, pkz, peln, q_con, ak, bk, & ks, gridstruct, flagstruct, neststruct, idiag, bd, domain, & init_step, i_pack, end_step, diss_est,time_total) @@ -244,7 +244,6 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, real, intent(inout):: vc(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz) real, intent(inout), dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz):: ua, va !The following 2 variables are for SA-3D-TKE (kyf) (modify for data structure) - logical, intent(IN) :: sa3dtke_dyco type(sa3dtke_type), intent(inout) :: sa3dtke_var real, intent(inout):: q_con(bd%isd:, bd%jsd:, 1:) @@ -301,8 +300,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, real :: split_timestep_bc !The following is for SA-3D-TKE - real:: tke(bd%isd:bd%ied,bd%jsd:bd%jed,npz) - integer :: ntke + integer :: sgs_tke integer :: is, ie, js, je integer :: isd, ied, jsd, jed @@ -782,14 +780,14 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, endif !The following is modified for SA-3D-TKE -!-- add sa3dtke_dyco and dku3d_h in parallel calculation +!-- add sa3dtke_var in parallel calculation ! replce dku3d_h by sa3dtke_var (kyf) (modify for data structure) call timing_on('d_sw') !$OMP parallel do default(none) shared(npz,flagstruct,nord_v,pfull,damp_vt,hydrostatic,last_step, & !$OMP is,ie,js,je,isd,ied,jsd,jed,omga,delp,gridstruct,npx,npy, & !$OMP ng,zh,vt,ptc,pt,u,v,w,uc,vc,ua,va,divgd,mfx,mfy,cx,cy, & -!$OMP sa3dtke_dyco, sa3dtke_var, & +!$OMP sa3dtke_var, & !$OMP crx,cry,xfx,yfx,q_con,zvir,sphum,nq,q,dt,bd,rdt,iep1,jep1, & !$OMP heat_source,diss_est,ptop,first_call) & !$OMP private(nord_k, nord_w, nord_t, damp_w, damp_t, d2_divg, & @@ -903,8 +901,8 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, u(isd,jsd,k), v(isd,jsd,k), w(isd:,jsd:,k), uc(isd,jsd,k), & vc(isd,jsd,k), ua(isd,jsd,k), va(isd,jsd,k), divgd(isd,jsd,k), & mfx(is, js, k), mfy(is, js, k), cx(is, jsd,k), cy(isd,js, k), & -!The following 2 variables are for SA-3D-TKE (kyf) (modify for data structure) - sa3dtke_dyco, sa3dtke_var%dku3d_h(isd, jsd, k), & +!The following variable is for SA-3D-TKE (kyf) (modify for data structure) + sa3dtke_var%dku3d_h(isd, jsd, k), & crx(is, jsd,k), cry(isd,js, k), xfx(is, jsd,k), yfx(isd,js, k), & #ifdef USE_COND q_con(isd:,jsd:,k), z_rat(isd,jsd), & @@ -1392,21 +1390,21 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, enddo ! time split loop first_call=.false. !----------------------------------------------------- -!The following is for SA-3D-TKE (kyf) (modify for data structure) +!The following is for SA-3D-TKE (kyf) (modify for data structure)(update mpi use) !--calculating shear deformation and TKE transport for 3d TKE scheme - if(sa3dtke_dyco) then + if(flagstruct%sa3dtke_dyco) then if ( end_step ) then - ntke = get_tracer_index(MODEL_ATMOS, 'sgs_tke') - call mpp_update_domains(q(:,:,:,ntke), domain) - call mpp_update_domains(sa3dtke_var%dku3d_e(:,:,:), domain) ! update dku3d_e at halo + sgs_tke = get_tracer_index(MODEL_ATMOS, 'sgs_tke') + call mpp_update_domains(q(:,:,:,sgs_tke), domain, complete=.false.) + call mpp_update_domains(sa3dtke_var%dku3d_e(:,:,:), domain, complete=.true.) ! update dku3d_e at halo call diff3d(npx, npy, npz, nq, ua, va, w, & q, sa3dtke_var%deform_1, sa3dtke_var%deform_2, & sa3dtke_var%deform_3, sa3dtke_var%dku3d_e, & zh, dp_ref, gridstruct, bd) - call mpp_update_domains(sa3dtke_var%deform_1, domain) - call mpp_update_domains(sa3dtke_var%deform_2, domain) - call mpp_update_domains(sa3dtke_var%deform_3, domain) + call mpp_update_domains(sa3dtke_var%deform_1, domain, complete=.false.) + call mpp_update_domains(sa3dtke_var%deform_2, domain, complete=.false.) + call mpp_update_domains(sa3dtke_var%deform_3, domain, complete=.true.) endif endif @@ -2199,7 +2197,7 @@ subroutine diff3d(npx, npy, npz, nq, ua, va, w, q, & real:: dkdx(bd%isd:bd%ied,bd%jsd:bd%jed,npz) real:: dkdy(bd%isd:bd%ied,bd%jsd:bd%jed,npz) - integer :: ntke + integer :: sgs_tke real :: tke_1(bd%isd:bd%ied+2, bd%jsd:bd%jed+2) real :: dedy_1(bd%isd:bd%ied+1, bd%jsd:bd%jed+1, npz) real :: dedy_2(bd%isd:bd%ied, bd%jsd:bd%jed, npz) @@ -2342,18 +2340,18 @@ subroutine diff3d(npx, npy, npz, nq, ua, va, w, q, & ! where e is tke !=========================================================== - ntke = get_tracer_index(MODEL_ATMOS, 'sgs_tke') + sgs_tke = get_tracer_index(MODEL_ATMOS, 'sgs_tke') ! KGao: make the 2d temporay arrays private - as suggested by Lucas !$OMP parallel do default(none) shared(npz,is,ie,js,je,ua,va,w,q,dx,dy,rarea, & -!$OMP dku3d_e,dkdx,dkdy,dedx_1,dedy_1,ntke,dedy_2,dedx_2) & +!$OMP dku3d_e,dkdx,dkdy,dedx_1,dedy_1,sgs_tke,dedy_2,dedx_2) & !$OMP private(ut,vt,tke_1) do k=1,npz !------------------------------------- do j=js,je+2 do i=is,ie+2 - tke_1(i,j)=q(i,j,k,ntke) + tke_1(i,j)=q(i,j,k,sgs_tke) enddo enddo do j=js,je+2 @@ -2422,7 +2420,7 @@ subroutine diff3d(npx, npy, npz, nq, ua, va, w, q, & ! get deform_3 based on all terms !$OMP parallel do default(none) shared(npz,is,ie,js,je,deform_3, & -!$OMP dku3d_e,dkdx,dkdy,dedx_1,dedy_1,ntke,dedx_2,dedy_2) +!$OMP dku3d_e,dkdx,dkdy,dedx_1,dedy_1,sgs_tke,dedx_2,dedy_2) do k=1,npz do j=js,je do i=is,ie diff --git a/model/fv_dynamics.F90 b/model/fv_dynamics.F90 index fc98e9b6c..063589282 100644 --- a/model/fv_dynamics.F90 +++ b/model/fv_dynamics.F90 @@ -189,8 +189,6 @@ module fv_dynamics_mod subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, & reproduce_sum, kappa, cp_air, zvir, ptop, ks, ncnst, n_split, & q_split, u, v, w, delz, hydrostatic, & -!The following variable is for SA-3D-TKE - sa3dtke_dyco, & pt, delp, q, & ps, pe, pk, peln, pkz, phis, q_con, omga, ua, va, uc, vc, & !The following variable is for SA-3D-TKE (kyf) (modify for data structure) @@ -262,7 +260,6 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, real, intent(inout), dimension(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz):: ua, va !The following variable is for SA-3D-TKE (kyf) (modify for data structure) - logical, intent(in) :: sa3dtke_dyco real, intent(in), dimension(npz+1):: ak, bk ! For SA-3D-TKE (kyf) (modify for data structure) @@ -677,8 +674,8 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill, grav, hydrostatic, & u, v, w, delz, pt, q, delp, pe, pk, phis, ws, omga, ptop, pfull, ua, va, & uc, vc, & -!The following 2 variables are for SA-3D-TKE (kyf) (modify for data structure) - sa3dtke_dyco, sa3dtke_var, & +!The following variable is for SA-3D-TKE (kyf) (modify for data structure) + sa3dtke_var, & mfx, mfy, cx, cy, pkz, peln, q_con, ak, bk, ks, & gridstruct, flagstruct, neststruct, idiag, bd, & domain, n_map==1, i_pack, GFDL_interstitial%last_step, diss_est,time_total) diff --git a/model/sw_core.F90 b/model/sw_core.F90 index 8e3bf0c85..866068a36 100644 --- a/model/sw_core.F90 +++ b/model/sw_core.F90 @@ -523,7 +523,7 @@ end subroutine c_sw subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, & ua, va, divg_d, xflux, yflux, cx, cy, & !The following 2 variables are for SA-3D-TKE (kyf) (modify for data structure) - sa3dtke_dyco, dku3d_h, & + dku3d_h, & crx_adv, cry_adv, xfx_adv, yfx_adv, q_con, z_rat, kgb, heat_source,diss_est, & zvir, sphum, nq, q, k, km, inline_q, & dt, hord_tr, hord_mt, hord_vt, hord_tm, hord_dp, nord, & @@ -541,7 +541,6 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, & real, intent(in):: damp_v, damp_w, damp_t, kgb type(fv_grid_bounds_type), intent(IN) :: bd !The following 2 variables are for SA-3D-TKE (kyf) (modify for data structure) - logical, intent(IN):: sa3dtke_dyco !type(sa3dtke_type), intent(in) :: sa3dtke_var real, intent(IN), dimension(bd%isd:bd%ied, bd%jsd:bd%jed):: dku3d_h @@ -1462,7 +1461,7 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, & do j=js,je+1 do i=is,ie+1 !The following is for SA-3D-TKE - if(sa3dtke_dyco) then + if(flagstruct%sa3dtke_dyco) then damp2 = abs(dt)*dku3d_h(i,j) else damp2 = gridstruct%da_min_c*max(d2_bg, min(0.20, dddmp*vort(i,j))) ! del-2 From ec21a2a5e0e70d167041f75f1470c1b223aaf063 Mon Sep 17 00:00:00 2001 From: "Bin.Liu" Date: Fri, 27 Jun 2025 01:21:49 +0000 Subject: [PATCH 23/29] Change dku3d_h to be an optional input variable for d_sw subroutine. And it will be passed in only if the sa3dtke_dyco is ".true.". This fixes debug build and run issues when sa3dtke_dyco is ".false." (default). From @samuelkyfung and @BinLiu-NOAA. --- model/dyn_core.F90 | 44 +++++++++++++++++++++++++++++++------------- model/sw_core.F90 | 12 +++++------- 2 files changed, 36 insertions(+), 20 deletions(-) diff --git a/model/dyn_core.F90 b/model/dyn_core.F90 index 6eb8c6d1b..d1d60d106 100644 --- a/model/dyn_core.F90 +++ b/model/dyn_core.F90 @@ -897,22 +897,40 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, enddo enddo endif - call d_sw(vt(isd,jsd,k), delp(isd,jsd,k), ptc(isd,jsd,k), pt(isd,jsd,k), & - u(isd,jsd,k), v(isd,jsd,k), w(isd:,jsd:,k), uc(isd,jsd,k), & - vc(isd,jsd,k), ua(isd,jsd,k), va(isd,jsd,k), divgd(isd,jsd,k), & - mfx(is, js, k), mfy(is, js, k), cx(is, jsd,k), cy(isd,js, k), & -!The following variable is for SA-3D-TKE (kyf) (modify for data structure) - sa3dtke_var%dku3d_h(isd, jsd, k), & - crx(is, jsd,k), cry(isd,js, k), xfx(is, jsd,k), yfx(isd,js, k), & + + if(flagstruct%sa3dtke_dyco) then + call d_sw(vt(isd,jsd,k), delp(isd,jsd,k), ptc(isd,jsd,k), pt(isd,jsd,k), & + u(isd,jsd,k), v(isd,jsd,k), w(isd:,jsd:,k), uc(isd,jsd,k), & + vc(isd,jsd,k), ua(isd,jsd,k), va(isd,jsd,k), divgd(isd,jsd,k), & + mfx(is, js, k), mfy(is, js, k), cx(is, jsd,k), cy(isd,js, k), & + crx(is, jsd,k), cry(isd,js, k), xfx(is, jsd,k), yfx(isd,js, k), & +#ifdef USE_COND + q_con(isd:,jsd:,k), z_rat(isd,jsd), & +#else + q_con(isd:,jsd:,1), z_rat(isd,jsd), & +#endif + kgb, heat_s, diss_e,zvir, sphum, nq, q, k, npz, flagstruct%inline_q, dt, & + flagstruct%hord_tr, hord_m, hord_v, hord_t, hord_p, & + nord_k, nord_v(k), nord_w, nord_t, flagstruct%dddmp, d2_divg, flagstruct%d4_bg, & + damp_vt(k), damp_w, damp_t, d_con_k, hydrostatic, gridstruct, flagstruct, bd, & +!The following optional variable is for SA-3D-TKE (kyf) + sa3dtke_var%dku3d_h(isd, jsd, k) ) + else + call d_sw(vt(isd,jsd,k), delp(isd,jsd,k), ptc(isd,jsd,k), pt(isd,jsd,k), & + u(isd,jsd,k), v(isd,jsd,k), w(isd:,jsd:,k), uc(isd,jsd,k), & + vc(isd,jsd,k), ua(isd,jsd,k), va(isd,jsd,k), divgd(isd,jsd,k), & + mfx(is, js, k), mfy(is, js, k), cx(is, jsd,k), cy(isd,js, k), & + crx(is, jsd,k), cry(isd,js, k), xfx(is, jsd,k), yfx(isd,js, k), & #ifdef USE_COND - q_con(isd:,jsd:,k), z_rat(isd,jsd), & + q_con(isd:,jsd:,k), z_rat(isd,jsd), & #else - q_con(isd:,jsd:,1), z_rat(isd,jsd), & + q_con(isd:,jsd:,1), z_rat(isd,jsd), & #endif - kgb, heat_s, diss_e,zvir, sphum, nq, q, k, npz, flagstruct%inline_q, dt, & - flagstruct%hord_tr, hord_m, hord_v, hord_t, hord_p, & - nord_k, nord_v(k), nord_w, nord_t, flagstruct%dddmp, d2_divg, flagstruct%d4_bg, & - damp_vt(k), damp_w, damp_t, d_con_k, hydrostatic, gridstruct, flagstruct, bd) + kgb, heat_s, diss_e,zvir, sphum, nq, q, k, npz, flagstruct%inline_q, dt, & + flagstruct%hord_tr, hord_m, hord_v, hord_t, hord_p, & + nord_k, nord_v(k), nord_w, nord_t, flagstruct%dddmp, d2_divg, flagstruct%d4_bg, & + damp_vt(k), damp_w, damp_t, d_con_k, hydrostatic, gridstruct, flagstruct, bd) + endif if((.not.flagstruct%use_old_omega) .and. last_step ) then ! Average horizontal "convergence" to cell center diff --git a/model/sw_core.F90 b/model/sw_core.F90 index 866068a36..54d7e7bd7 100644 --- a/model/sw_core.F90 +++ b/model/sw_core.F90 @@ -50,7 +50,6 @@ module sw_core_mod deln_flux_explm, deln_flux_explm_udvd use fv_mp_mod, only: is_master, fill_corners, XDir, YDir use fv_arrays_mod, only: fv_grid_type, fv_grid_bounds_type, fv_flags_type - !use fv_arrays_mod, only: sa3dtke_type ! for SA-3D-TKE (kyf) (modify for data structure) use a2b_edge_mod, only: a2b_ord4 #ifdef SW_DYNAMICS @@ -522,13 +521,13 @@ end subroutine c_sw !! and other prognostic varaiables. subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, & ua, va, divg_d, xflux, yflux, cx, cy, & -!The following 2 variables are for SA-3D-TKE (kyf) (modify for data structure) - dku3d_h, & crx_adv, cry_adv, xfx_adv, yfx_adv, q_con, z_rat, kgb, heat_source,diss_est, & zvir, sphum, nq, q, k, km, inline_q, & dt, hord_tr, hord_mt, hord_vt, hord_tm, hord_dp, nord, & nord_v, nord_w, nord_t, dddmp, d2_bg, d4_bg, damp_v, damp_w, & - damp_t, d_con, hydrostatic, gridstruct, flagstruct, bd) + damp_t, d_con, hydrostatic, gridstruct, flagstruct, bd, & +!The following optional variable is for SA-3D-TKE (kyf) + dku3d_h) integer, intent(IN):: hord_tr, hord_mt, hord_vt, hord_tm, hord_dp integer, intent(IN):: nord !< nord=1 divergence damping; (del-4) or 3 (del-8) @@ -540,9 +539,6 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, & real , intent(IN):: zvir real, intent(in):: damp_v, damp_w, damp_t, kgb type(fv_grid_bounds_type), intent(IN) :: bd -!The following 2 variables are for SA-3D-TKE (kyf) (modify for data structure) - !type(sa3dtke_type), intent(in) :: sa3dtke_var - real, intent(IN), dimension(bd%isd:bd%ied, bd%jsd:bd%jed):: dku3d_h real, intent(inout):: divg_d(bd%isd:bd%ied+1,bd%jsd:bd%jed+1) !< divergence real, intent(IN), dimension(bd%isd:bd%ied, bd%jsd:bd%jed):: z_rat @@ -566,6 +562,8 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, & real, intent(OUT), dimension(bd%isd:bd%ied,bd%js:bd%je+1):: cry_adv, yfx_adv type(fv_grid_type), intent(IN), target :: gridstruct type(fv_flags_type), intent(IN), target :: flagstruct +!The following optional variable is for SA-3D-TKE (kyf) + real, intent(IN), optional, dimension(bd%isd:bd%ied,bd%jsd:bd%jed):: dku3d_h ! Local: logical:: sw_corner, se_corner, ne_corner, nw_corner real :: ut(bd%isd:bd%ied+1,bd%jsd:bd%jed) From d7b9939b6812f7d8b06e87807b4968cec301df0e Mon Sep 17 00:00:00 2001 From: Bin Liu Date: Thu, 7 Aug 2025 08:48:36 +0000 Subject: [PATCH 24/29] From @JongilHan66: For the 3DTKE option, update damp2 (in sw_core.F90) by using the maximum of Smagorinsky divergence damping and dt*Kh (from 3dtke). Basically, damp2=max([Smagorinsky damping],dt*Kh). Note: Based on the testing with HAFS for some 2023-2024 hurricanes, using max([Smagorinsky damping],dt*Kh) instead of sum([Smagorinsky damping],dt*Kh) for damp2 provides overall better intensity vmax biases with slightly better intensity forecast skills at longer forecast lead times of days 4-5. --- model/sw_core.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/model/sw_core.F90 b/model/sw_core.F90 index eef80f387..53b1d1c02 100644 --- a/model/sw_core.F90 +++ b/model/sw_core.F90 @@ -1461,7 +1461,7 @@ subroutine d_sw(delpc, delp, ptc, pt, u, v, w, uc,vc, & damp2 = gridstruct%da_min_c*max(d2_bg, min(0.20, dddmp*vort(i,j))) ! del-2 !The following is for SA-3D-TKE if(flagstruct%sa3dtke_dyco) then - damp2 = damp2 + abs(dt)*dku3d_h(i,j) + damp2 = max(damp2, abs(dt)*dku3d_h(i,j)) endif vort(i,j) = damp2*delpc(i,j) + dd8*divg_d(i,j) From 9660070f42bb945dd5d9a620fb3283a87975eed2 Mon Sep 17 00:00:00 2001 From: XuLu-NOAA Date: Wed, 10 Dec 2025 17:27:58 +0000 Subject: [PATCH 25/29] Using the FV3 read_cubed_sphere_inc function to read in the regional IAU increments. --- tools/cubed_sphere_inc_mod.F90 | 84 +++++++++++++++++++++++----------- tools/fv_iau_mod.F90 | 50 +++++--------------- tools/fv_treat_da_inc.F90 | 2 +- 3 files changed, 70 insertions(+), 66 deletions(-) diff --git a/tools/cubed_sphere_inc_mod.F90 b/tools/cubed_sphere_inc_mod.F90 index 2b169b6e7..aa6c1a0fb 100644 --- a/tools/cubed_sphere_inc_mod.F90 +++ b/tools/cubed_sphere_inc_mod.F90 @@ -4,8 +4,9 @@ module cubed_sphere_inc_mod use field_manager_mod, only: MODEL_ATMOS use fv_arrays_mod, only: fv_atmos_type use fms2_io_mod, only: open_file, close_file, read_data, variable_exists, & - FmsNetcdfDomainFile_t, register_axis - + FmsNetcdfDomainFile_t, register_axis, FmsNetcdfFile_t + use mpp_domains_mod, only: mpp_get_compute_domain + implicit none type increment_data_type real, allocatable :: ua_inc(:,:,:) @@ -22,44 +23,73 @@ module cubed_sphere_inc_mod !---------------------------------------------------------------------------------------- - subroutine read_cubed_sphere_inc(fname, increment_data, Atm) + subroutine read_cubed_sphere_inc(fname, increment_data, Atm, IAU_regional) character(*), intent(in) :: fname type(increment_data_type), intent(inout) :: increment_data type(fv_atmos_type), intent(in) :: Atm + logical, intent(in) :: IAU_regional type(FmsNetcdfDomainFile_t) :: fileobj + type(FmsNetcdfFile_t) :: fileobj_regional integer :: itracer, ntracers character(len=64) :: tracer_name - + integer :: is, ie, js, je + real:: tmp (Atm%npx-1,Atm%npy-1,Atm%npz) ! Get various dimensions call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers) + call mpp_get_compute_domain(Atm%domain, is, ie, js, je) ! Open file - if ( open_file(fileobj, trim(fname), 'read', Atm%domain) ) then - ! Register axes - call register_axis(fileobj, 'xaxis_1', 'x') - call register_axis(fileobj, 'yaxis_1', 'y') - - ! Read increments - call read_data(fileobj, 'u_inc', increment_data%ua_inc) - call read_data(fileobj, 'v_inc', increment_data%va_inc) - call read_data(fileobj, 'T_inc', increment_data%temp_inc) - call read_data(fileobj, 'delp_inc', increment_data%delp_inc) - if ( .not. Atm%flagstruct%hydrostatic ) then - call read_data(fileobj, 'delz_inc', increment_data%delz_inc) - end if + if (IAU_regional) then + if ( open_file(fileobj_regional, trim(fname), 'read') ) then + ! Read increments + call read_data(fileobj_regional, 'u_inc', tmp) + increment_data%ua_inc=tmp(is:ie,js:je,:) + call read_data(fileobj_regional, 'v_inc', tmp) + increment_data%va_inc=tmp(is:ie,js:je,:) + call read_data(fileobj_regional, 'T_inc', tmp) + increment_data%temp_inc=tmp(is:ie,js:je,:) + call read_data(fileobj_regional, 'delp_inc', tmp) + increment_data%delp_inc=tmp(is:ie,js:je,:) + if ( .not. Atm%flagstruct%hydrostatic ) then + call read_data(fileobj_regional, 'delz_inc', tmp) + increment_data%delz_inc=tmp(is:ie,js:je,:) + end if + + ! Read tracer increments + do itracer = 1,ntracers + call get_tracer_names(MODEL_ATMOS, itracer, tracer_name) + if ( variable_exists(fileobj_regional, trim(tracer_name)//'_inc') ) then + call read_data(fileobj_regional, trim(tracer_name)//'_inc', tmp) + increment_data%tracer_inc(:,:,:,itracer)=tmp(is:ie,js:je,:) + end if + end do + call close_file(fileobj_regional) + end if + else + if ( open_file(fileobj, trim(fname), 'read', Atm%domain) ) then + ! Register axes + call register_axis(fileobj, 'xaxis_1', 'x') + call register_axis(fileobj, 'yaxis_1', 'y') + ! Read increments + call read_data(fileobj, 'u_inc', increment_data%ua_inc) + call read_data(fileobj, 'v_inc', increment_data%va_inc) + call read_data(fileobj, 'T_inc', increment_data%temp_inc) + call read_data(fileobj, 'delp_inc', increment_data%delp_inc) + if ( .not. Atm%flagstruct%hydrostatic ) then + call read_data(fileobj, 'delz_inc', increment_data%delz_inc) + end if - ! Read tracer increments - do itracer = 1,ntracers - call get_tracer_names(MODEL_ATMOS, itracer, tracer_name) - if ( variable_exists(fileobj, trim(tracer_name)//'_inc') ) then - call read_data(fileobj, trim(tracer_name)//'_inc', increment_data%tracer_inc(:,:,:,itracer)) - end if - end do - - call close_file(fileobj) + ! Read tracer increments + do itracer = 1,ntracers + call get_tracer_names(MODEL_ATMOS, itracer, tracer_name) + if ( variable_exists(fileobj, trim(tracer_name)//'_inc') ) then + call read_data(fileobj, trim(tracer_name)//'_inc', increment_data%tracer_inc(:,:,:,itracer)) + end if + end do + call close_file(fileobj) + end if end if - end subroutine read_cubed_sphere_inc !---------------------------------------------------------------------------------------- diff --git a/tools/fv_iau_mod.F90 b/tools/fv_iau_mod.F90 index 09da2af1a..2b464f153 100644 --- a/tools/fv_iau_mod.F90 +++ b/tools/fv_iau_mod.F90 @@ -280,7 +280,7 @@ subroutine IAU_initialize (IPD_Control, IAU_Data, Init_parm, Atm) iau_state%wt_normfact = (2*nstep+1)/normfact endif if ( Atm%flagstruct%increment_file_on_native_grid ) then - call read_cubed_sphere_inc('INPUT/'//trim(IPD_Control%iau_inc_files(1)), iau_state%inc1, Atm) + call read_cubed_sphere_inc('INPUT/'//trim(IPD_Control%iau_inc_files(1)), iau_state%inc1, Atm, IPD_Control%iau_regional) else call read_iau_forcing(IPD_Control,iau_state%inc1,'INPUT/'//trim(IPD_Control%iau_inc_files(1))) endif @@ -296,7 +296,7 @@ subroutine IAU_initialize (IPD_Control, IAU_Data, Init_parm, Atm) allocate (iau_state%inc2%tracer_inc(is:ie, js:je, km,ntracers)) iau_state%hr2=IPD_Control%iaufhrs(2) if ( Atm%flagstruct%increment_file_on_native_grid ) then - call read_cubed_sphere_inc('INPUT/'//trim(IPD_Control%iau_inc_files(2)), iau_state%inc2, Atm) + call read_cubed_sphere_inc('INPUT/'//trim(IPD_Control%iau_inc_files(2)), iau_state%inc2, Atm, IPD_Control%iau_regional) else call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(2))) endif @@ -386,7 +386,7 @@ subroutine getiauforcing(IPD_Control,IAU_Data,Atm) iau_state%inc1=iau_state%inc2 if (is_master()) print *,'reading next increment file',trim(IPD_Control%iau_inc_files(itnext)) if ( Atm%flagstruct%increment_file_on_native_grid ) then - call read_cubed_sphere_inc('INPUT/'//trim(IPD_Control%iau_inc_files(itnext)), iau_state%inc2, Atm) + call read_cubed_sphere_inc('INPUT/'//trim(IPD_Control%iau_inc_files(itnext)), iau_state%inc2, Atm, IPD_Control%iau_regional) else call read_iau_forcing(IPD_Control,iau_state%inc2,'INPUT/'//trim(IPD_Control%iau_inc_files(itnext))) endif @@ -407,6 +407,7 @@ subroutine updateiauforcing(IPD_Control,IAU_Data,wt) ! if (is_master()) print *,'in updateiauforcing',nfiles,IPD_Control%iaufhrs(1:nfiles) delt = (iau_state%hr2-(IPD_Control%fhour))/(IAU_state%hr2-IAU_state%hr1) + if (is_master()) print *,'XL',maxval(IAU_state%inc2%ua_inc(is:ie,js:je,:)),minval(IAU_state%inc2%ua_inc(is:ie,js:je,:)),delt,IPD_Control%iau_regional do j = js,je do i = is,ie do k = 1,npz @@ -499,25 +500,14 @@ subroutine read_iau_forcing(IPD_Control,increments,fname) endif ! read in 1 time level - if (IPD_Control%iau_regional) then - call interp_inc2('T_inc',increments%temp_inc) - call interp_inc2('delp_inc',increments%delp_inc) - call interp_inc2('delz_inc',increments%delz_inc) - call interp_inc2('u_inc',increments%ua_inc) - call interp_inc2('v_inc',increments%va_inc) - do l=1,ntracers - call interp_inc2(trim(tracer_names(l))//'_inc',increments%tracer_inc(:,:,:,l)) - enddo - else - call interp_inc('T_inc',increments%temp_inc(:,:,:),jbeg,jend) - call interp_inc('delp_inc',increments%delp_inc(:,:,:),jbeg,jend) - call interp_inc('delz_inc',increments%delz_inc(:,:,:),jbeg,jend) - call interp_inc('u_inc',increments%ua_inc(:,:,:),jbeg,jend) ! can these be treated as scalars? - call interp_inc('v_inc',increments%va_inc(:,:,:),jbeg,jend) - do l=1,ntracers - call interp_inc(trim(tracer_names(l))//'_inc',increments%tracer_inc(:,:,:,l),jbeg,jend) - enddo - end if + call interp_inc('T_inc',increments%temp_inc(:,:,:),jbeg,jend) + call interp_inc('delp_inc',increments%delp_inc(:,:,:),jbeg,jend) + call interp_inc('delz_inc',increments%delz_inc(:,:,:),jbeg,jend) + call interp_inc('u_inc',increments%ua_inc(:,:,:),jbeg,jend) ! can these be treated as scalars? + call interp_inc('v_inc',increments%va_inc(:,:,:),jbeg,jend) + do l=1,ntracers + call interp_inc(trim(tracer_names(l))//'_inc',increments%tracer_inc(:,:,:,l),jbeg,jend) + enddo call close_ncfile(ncid) deallocate (wk3) @@ -551,22 +541,6 @@ subroutine interp_inc(field_name,var,jbeg,jend) enddo end subroutine interp_inc -subroutine interp_inc2(field_name,var) -! interpolate increment from GSI gaussian grid to cubed sphere -! everying is on the A-grid, earth relative - character(len=*), intent(in) :: field_name - real, dimension(is:ie,js:je,1:km), intent(inout) :: var - integer:: i1, i2, j1, k,j,i,ierr - call check_var_exists(ncid, field_name, ierr) - if (ierr == 0) then - call get_var3_r4( ncid, field_name, 1,im, js,je, 1,km, wk3 ) - else - if (is_master()) print *,'warning: no increment for ',trim(field_name),' found, assuming zero' - wk3 = 0. - end if - var(is:ie,js:je,1:km)=wk3(is:ie,:,:) -end subroutine interp_inc2 - end module fv_iau_mod diff --git a/tools/fv_treat_da_inc.F90 b/tools/fv_treat_da_inc.F90 index 7e985fb13..6d749bf15 100644 --- a/tools/fv_treat_da_inc.F90 +++ b/tools/fv_treat_da_inc.F90 @@ -522,7 +522,7 @@ subroutine read_da_inc_cubed_sphere(Atm, fv_domain, bd, npz_in, nq, & ! Read increments fname = trim(fname_prefix) // '.nc' - call read_cubed_sphere_inc(fname, increment_data, Atm) + call read_cubed_sphere_inc(fname, increment_data, Atm, .False.) ! Wind increments ! --------------- From 147de571066703078fe3a6972c55d9a22c73f45f Mon Sep 17 00:00:00 2001 From: XuLu-NOAA Date: Wed, 10 Dec 2025 21:44:18 +0000 Subject: [PATCH 26/29] Change Regional 3DIAU from even weighting to linear decaying weighting --- tools/fv_iau_mod.F90 | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/tools/fv_iau_mod.F90 b/tools/fv_iau_mod.F90 index 2b464f153..4dbd37bdf 100644 --- a/tools/fv_iau_mod.F90 +++ b/tools/fv_iau_mod.F90 @@ -259,6 +259,11 @@ subroutine IAU_initialize (IPD_Control, IAU_Data, Init_parm, Atm) iau_state%wt = 1.0 ! IAU increment filter weights (default 1.0) iau_state%wt=iau_state%wt*IPD_Control%iau_inc_scale !Increase the weight iau_state%wt_normfact = 1.0 + if (IPD_Control%iau_regional) then + dtp=IPD_control%dtp + nstep = 0.5*IPD_Control%iau_delthrs*3600./dtp + iau_state%wt=(1.0/real(nstep))*(1.0-dtp/real(nstep))*IPD_Control%iau_inc_scale + end if if (IPD_Control%iau_filter_increments) then ! compute increment filter weights, sum to obtain normalization factor dtp=IPD_control%dtp @@ -359,8 +364,15 @@ subroutine getiauforcing(IPD_Control,IAU_Data,Atm) ! if (is_master()) print *,'no iau forcing',t1,IPD_Control%fhour,t2 IAU_Data%in_interval=.false. else + if (IPD_Control%iau_regional) then + dtp=IPD_control%dtp + nstep = 0.5*IPD_Control%iau_delthrs*3600./dtp + kstep = (IPD_Control%fhour-t1)*3600./dtp-nstep + iau_state%wt=(1.0/real(nstep))*(1.0-kstep/real(nstep))*IPD_Control%iau_inc_scale + call setiauforcing(IPD_Control,IAU_Data,iau_state%wt) + end if if (IPD_Control%iau_filter_increments) call setiauforcing(IPD_Control,IAU_Data,iau_state%wt) - if (is_master()) print *,'apply iau forcing t1,t,t2,filter wt=',t1,IPD_Control%fhour,t2,iau_state%wt/iau_state%wt_normfact,iau_state%wt + if (is_master()) print *,'apply iau forcing t1,t,t2,filter wt=',t1,IPD_Control%fhour,t2,iau_state%wt/iau_state%wt_normfact,iau_state%wt,kstep,nstep IAU_Data%in_interval=.true. endif return From 9fc595be672f2d829fa99fdf5cce2812362e1d9e Mon Sep 17 00:00:00 2001 From: XuLu-NOAA Date: Fri, 12 Dec 2025 03:17:07 +0000 Subject: [PATCH 27/29] Revert the linear decaying weighting function in 3DIAU. Switch to the Lanczos-like temporal filter --- tools/fv_iau_mod.F90 | 19 ++++--------------- 1 file changed, 4 insertions(+), 15 deletions(-) diff --git a/tools/fv_iau_mod.F90 b/tools/fv_iau_mod.F90 index 4dbd37bdf..0ce8277b1 100644 --- a/tools/fv_iau_mod.F90 +++ b/tools/fv_iau_mod.F90 @@ -259,15 +259,10 @@ subroutine IAU_initialize (IPD_Control, IAU_Data, Init_parm, Atm) iau_state%wt = 1.0 ! IAU increment filter weights (default 1.0) iau_state%wt=iau_state%wt*IPD_Control%iau_inc_scale !Increase the weight iau_state%wt_normfact = 1.0 - if (IPD_Control%iau_regional) then - dtp=IPD_control%dtp - nstep = 0.5*IPD_Control%iau_delthrs*3600./dtp - iau_state%wt=(1.0/real(nstep))*(1.0-dtp/real(nstep))*IPD_Control%iau_inc_scale - end if if (IPD_Control%iau_filter_increments) then ! compute increment filter weights, sum to obtain normalization factor dtp=IPD_control%dtp - nstep = 0.5*IPD_Control%iau_delthrs*3600/dtp + nstep = nint(0.5*IPD_Control%iau_delthrs*3600/dtp) ! compute normalization factor for filter weights normfact = 0. do k=1,2*nstep+1 @@ -284,6 +279,7 @@ subroutine IAU_initialize (IPD_Control, IAU_Data, Init_parm, Atm) enddo iau_state%wt_normfact = (2*nstep+1)/normfact endif + if (IPD_Control%iau_regional) iau_state%wt_normfact=iau_state%wt_normfact*IPD_Control%iau_inc_scale if ( Atm%flagstruct%increment_file_on_native_grid ) then call read_cubed_sphere_inc('INPUT/'//trim(IPD_Control%iau_inc_files(1)), iau_state%inc1, Atm, IPD_Control%iau_regional) else @@ -339,9 +335,9 @@ subroutine getiauforcing(IPD_Control,IAU_Data,Atm) ! in window kstep=-nstep,nstep (2*nstep+1 total) ! time step IPD_control%dtp dtp=IPD_control%dtp - nstep = 0.5*IPD_Control%iau_delthrs*3600/dtp + nstep = nint(0.5*IPD_Control%iau_delthrs*3600/dtp) ! compute normalized filter weight - kstep = ((IPD_Control%fhour-t1) - 0.5*IPD_Control%iau_delthrs)*3600./dtp + kstep = nint(((IPD_Control%fhour-t1) - 0.5*IPD_Control%iau_delthrs)*3600./dtp) if (IPD_Control%fhour >= t1 .and. IPD_Control%fhour < t2) then sx = acos(-1.)*kstep/nstep wx = acos(-1.)*kstep/(nstep+1) @@ -364,13 +360,6 @@ subroutine getiauforcing(IPD_Control,IAU_Data,Atm) ! if (is_master()) print *,'no iau forcing',t1,IPD_Control%fhour,t2 IAU_Data%in_interval=.false. else - if (IPD_Control%iau_regional) then - dtp=IPD_control%dtp - nstep = 0.5*IPD_Control%iau_delthrs*3600./dtp - kstep = (IPD_Control%fhour-t1)*3600./dtp-nstep - iau_state%wt=(1.0/real(nstep))*(1.0-kstep/real(nstep))*IPD_Control%iau_inc_scale - call setiauforcing(IPD_Control,IAU_Data,iau_state%wt) - end if if (IPD_Control%iau_filter_increments) call setiauforcing(IPD_Control,IAU_Data,iau_state%wt) if (is_master()) print *,'apply iau forcing t1,t,t2,filter wt=',t1,IPD_Control%fhour,t2,iau_state%wt/iau_state%wt_normfact,iau_state%wt,kstep,nstep IAU_Data%in_interval=.true. From b40d4cd2613bfbc4f18b2038576c75ca175058c0 Mon Sep 17 00:00:00 2001 From: XuLu-NOAA Date: Mon, 29 Dec 2025 14:44:44 +0000 Subject: [PATCH 28/29] Adjust the regional 3DIAU interface with the global FMS approach with the help of Uriel Ramirez (uriel.ramirez@noaa.gov) --- tools/cubed_sphere_inc_mod.F90 | 85 ++++++++++++++-------------------- tools/fv_iau_mod.F90 | 6 +-- 2 files changed, 38 insertions(+), 53 deletions(-) diff --git a/tools/cubed_sphere_inc_mod.F90 b/tools/cubed_sphere_inc_mod.F90 index aa6c1a0fb..2e0dc8fa4 100644 --- a/tools/cubed_sphere_inc_mod.F90 +++ b/tools/cubed_sphere_inc_mod.F90 @@ -5,7 +5,7 @@ module cubed_sphere_inc_mod use fv_arrays_mod, only: fv_atmos_type use fms2_io_mod, only: open_file, close_file, read_data, variable_exists, & FmsNetcdfDomainFile_t, register_axis, FmsNetcdfFile_t - use mpp_domains_mod, only: mpp_get_compute_domain + use fv_mp_mod, only: is_master implicit none type increment_data_type @@ -33,62 +33,47 @@ subroutine read_cubed_sphere_inc(fname, increment_data, Atm, IAU_regional) type(FmsNetcdfFile_t) :: fileobj_regional integer :: itracer, ntracers character(len=64) :: tracer_name - integer :: is, ie, js, je - real:: tmp (Atm%npx-1,Atm%npy-1,Atm%npz) + character(len=:), allocatable :: fname_base + integer :: ipos ! Get various dimensions call get_number_tracers(MODEL_ATMOS, num_tracers=ntracers) - call mpp_get_compute_domain(Atm%domain, is, ie, js, je) ! Open file if (IAU_regional) then - if ( open_file(fileobj_regional, trim(fname), 'read') ) then - ! Read increments - call read_data(fileobj_regional, 'u_inc', tmp) - increment_data%ua_inc=tmp(is:ie,js:je,:) - call read_data(fileobj_regional, 'v_inc', tmp) - increment_data%va_inc=tmp(is:ie,js:je,:) - call read_data(fileobj_regional, 'T_inc', tmp) - increment_data%temp_inc=tmp(is:ie,js:je,:) - call read_data(fileobj_regional, 'delp_inc', tmp) - increment_data%delp_inc=tmp(is:ie,js:je,:) - if ( .not. Atm%flagstruct%hydrostatic ) then - call read_data(fileobj_regional, 'delz_inc', tmp) - increment_data%delz_inc=tmp(is:ie,js:je,:) - end if - - ! Read tracer increments - do itracer = 1,ntracers - call get_tracer_names(MODEL_ATMOS, itracer, tracer_name) - if ( variable_exists(fileobj_regional, trim(tracer_name)//'_inc') ) then - call read_data(fileobj_regional, trim(tracer_name)//'_inc', tmp) - increment_data%tracer_inc(:,:,:,itracer)=tmp(is:ie,js:je,:) - end if - end do - call close_file(fileobj_regional) - end if + ipos = index(fname, '.') + if (ipos > 0) then + fname_base = fname(1:ipos-1) + else + fname_base = trim(fname) + endif + fname_base=trim(fname_base)//'.nc' else - if ( open_file(fileobj, trim(fname), 'read', Atm%domain) ) then - ! Register axes - call register_axis(fileobj, 'xaxis_1', 'x') - call register_axis(fileobj, 'yaxis_1', 'y') - ! Read increments - call read_data(fileobj, 'u_inc', increment_data%ua_inc) - call read_data(fileobj, 'v_inc', increment_data%va_inc) - call read_data(fileobj, 'T_inc', increment_data%temp_inc) - call read_data(fileobj, 'delp_inc', increment_data%delp_inc) - if ( .not. Atm%flagstruct%hydrostatic ) then - call read_data(fileobj, 'delz_inc', increment_data%delz_inc) - end if + fname_base=trim(fname) + end if + if ( open_file(fileobj, trim(fname_base), 'read', Atm%domain) ) then + ! Register axes + call register_axis(fileobj, 'xaxis_1', 'x') + call register_axis(fileobj, 'yaxis_1', 'y') + ! Read increments + call read_data(fileobj, 'u_inc', increment_data%ua_inc) + call read_data(fileobj, 'v_inc', increment_data%va_inc) + call read_data(fileobj, 'T_inc', increment_data%temp_inc) + call read_data(fileobj, 'delp_inc', increment_data%delp_inc) + if ( .not. Atm%flagstruct%hydrostatic ) then + call read_data(fileobj, 'delz_inc', increment_data%delz_inc) + end if - ! Read tracer increments - do itracer = 1,ntracers - call get_tracer_names(MODEL_ATMOS, itracer, tracer_name) - if ( variable_exists(fileobj, trim(tracer_name)//'_inc') ) then - call read_data(fileobj, trim(tracer_name)//'_inc', increment_data%tracer_inc(:,:,:,itracer)) - end if - end do - call close_file(fileobj) - end if + ! Read tracer increments + do itracer = 1,ntracers + call get_tracer_names(MODEL_ATMOS, itracer, tracer_name) + if ( variable_exists(fileobj, trim(tracer_name)//'_inc') ) then + call read_data(fileobj, trim(tracer_name)//'_inc', increment_data%tracer_inc(:,:,:,itracer)) + else + if (is_master()) print *,'warning: no increment for ',trim(tracer_name),' found, assuming zero' + increment_data%tracer_inc(:,:,:,itracer) = 0.0 + end if + end do + call close_file(fileobj) end if end subroutine read_cubed_sphere_inc diff --git a/tools/fv_iau_mod.F90 b/tools/fv_iau_mod.F90 index 0ce8277b1..d80a0c68d 100644 --- a/tools/fv_iau_mod.F90 +++ b/tools/fv_iau_mod.F90 @@ -183,9 +183,9 @@ subroutine IAU_initialize (IPD_Control, IAU_Data, Init_parm, Atm) if( file_exists(fname) ) then call open_ncfile( fname, ncid ) ! open the file if (IPD_Control%iau_regional) then - call get_ncdim1( ncid, 'nx', im) - call get_ncdim1( ncid, 'ny', jm) - call get_ncdim1( ncid, 'nz', km) + call get_ncdim1( ncid, 'xaxis_1', im) + call get_ncdim1( ncid, 'yaxis_1', jm) + call get_ncdim1( ncid, 'zaxis_1', km) else call get_ncdim1( ncid, 'lon', im) call get_ncdim1( ncid, 'lat', jm) From 3f742540a1beb0018dd4a54a603c3831e475a481 Mon Sep 17 00:00:00 2001 From: XuLu-NOAA Date: Mon, 12 Jan 2026 18:37:10 +0000 Subject: [PATCH 29/29] Changing is_master to mpp_error for warning messages --- tools/cubed_sphere_inc_mod.F90 | 4 ++-- tools/fv_iau_mod.F90 | 3 +-- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/tools/cubed_sphere_inc_mod.F90 b/tools/cubed_sphere_inc_mod.F90 index 2e0dc8fa4..19f807427 100644 --- a/tools/cubed_sphere_inc_mod.F90 +++ b/tools/cubed_sphere_inc_mod.F90 @@ -5,7 +5,7 @@ module cubed_sphere_inc_mod use fv_arrays_mod, only: fv_atmos_type use fms2_io_mod, only: open_file, close_file, read_data, variable_exists, & FmsNetcdfDomainFile_t, register_axis, FmsNetcdfFile_t - use fv_mp_mod, only: is_master + use mpp_mod, only: mpp_error, NOTE implicit none type increment_data_type @@ -69,7 +69,7 @@ subroutine read_cubed_sphere_inc(fname, increment_data, Atm, IAU_regional) if ( variable_exists(fileobj, trim(tracer_name)//'_inc') ) then call read_data(fileobj, trim(tracer_name)//'_inc', increment_data%tracer_inc(:,:,:,itracer)) else - if (is_master()) print *,'warning: no increment for ',trim(tracer_name),' found, assuming zero' + call mpp_error(NOTE, 'No Increment for '//trim(tracer_name)//' found, assuming zero') increment_data%tracer_inc(:,:,:,itracer) = 0.0 end if end do diff --git a/tools/fv_iau_mod.F90 b/tools/fv_iau_mod.F90 index d80a0c68d..38aa3e141 100644 --- a/tools/fv_iau_mod.F90 +++ b/tools/fv_iau_mod.F90 @@ -361,7 +361,7 @@ subroutine getiauforcing(IPD_Control,IAU_Data,Atm) IAU_Data%in_interval=.false. else if (IPD_Control%iau_filter_increments) call setiauforcing(IPD_Control,IAU_Data,iau_state%wt) - if (is_master()) print *,'apply iau forcing t1,t,t2,filter wt=',t1,IPD_Control%fhour,t2,iau_state%wt/iau_state%wt_normfact,iau_state%wt,kstep,nstep + if (is_master()) print *,'apply iau forcing t1,t,t2,filter wt=',t1,IPD_Control%fhour,t2,iau_state%wt/iau_state%wt_normfact IAU_Data%in_interval=.true. endif return @@ -408,7 +408,6 @@ subroutine updateiauforcing(IPD_Control,IAU_Data,wt) ! if (is_master()) print *,'in updateiauforcing',nfiles,IPD_Control%iaufhrs(1:nfiles) delt = (iau_state%hr2-(IPD_Control%fhour))/(IAU_state%hr2-IAU_state%hr1) - if (is_master()) print *,'XL',maxval(IAU_state%inc2%ua_inc(is:ie,js:je,:)),minval(IAU_state%inc2%ua_inc(is:ie,js:je,:)),delt,IPD_Control%iau_regional do j = js,je do i = is,ie do k = 1,npz