diff --git a/moving_nest/fv_moving_nest.F90 b/moving_nest/fv_moving_nest.F90 index 6ef5ab384..1c1288ca0 100644 --- a/moving_nest/fv_moving_nest.F90 +++ b/moving_nest/fv_moving_nest.F90 @@ -87,7 +87,7 @@ module fv_moving_nest_mod use fv_nwp_nudge_mod, only: do_adiabatic_init use init_hydro_mod, only: p_var use tracer_manager_mod, only: get_tracer_index, get_tracer_names - use fv_moving_nest_types_mod, only: fv_moving_nest_prog_type, fv_moving_nest_physics_type, Moving_nest + use fv_moving_nest_types_mod, only: fv_moving_nest_prog_type, fv_moving_nest_physics_type, Moving_nest, mn_land_mask_grids, mn_fix_grids, alloc_set_facwf use fv_moving_nest_utils_mod, only: alloc_halo_buffer, load_nest_latlons_from_nc, grid_geometry, output_grid_to_nc use fv_moving_nest_utils_mod, only: fill_nest_from_buffer, fill_nest_from_buffer_cell_center, fill_nest_from_buffer_nearest_neighbor use fv_moving_nest_utils_mod, only: fill_nest_halos_from_parent, fill_grid_from_supergrid, fill_weight_grid @@ -130,11 +130,13 @@ module fv_moving_nest_mod !! Step 6 interface mn_var_shift_data module procedure mn_var_shift_data_r4_2d - module procedure mn_var_shift_data_r4_3d + module procedure mn_var_shift_data_r4_3d_highz + module procedure mn_var_shift_data_r4_3d_lowhighz module procedure mn_var_shift_data_r4_4d module procedure mn_var_shift_data_r8_2d - module procedure mn_var_shift_data_r8_3d + module procedure mn_var_shift_data_r8_3d_highz + module procedure mn_var_shift_data_r8_3d_lowhighz module procedure mn_var_shift_data_r8_4d end interface mn_var_shift_data @@ -614,10 +616,10 @@ subroutine mn_latlon_load_parent(surface_dir, Atm, n, parent_tile, delta_i_c, de parent_geo%nxp = Atm(1)%npx parent_geo%nyp = Atm(1)%npy - + parent_geo%nx = parent_geo%nxp - 1 parent_geo%ny = parent_geo%nyp - 1 - + call mn_static_filename(surface_dir, parent_tile, 'grid', 1, grid_filename) call load_nest_latlons_from_nc(grid_filename, parent_geo%nxp, parent_geo%nyp, 1, pelist, & parent_geo, p_istart_fine, p_iend_fine, p_jstart_fine, p_jend_fine) @@ -803,6 +805,7 @@ subroutine mn_latlon_read_hires_parent(npx, npy, refine, pelist, fp_super_tile_g call load_nest_latlons_from_nc(trim(grid_filename), npx, npy, refine, pelist, & fp_super_tile_geo, fp_super_istart_fine, fp_super_iend_fine, fp_super_jstart_fine, fp_super_jend_fine) + end subroutine mn_latlon_read_hires_parent !>@brief The subroutine 'mn_orog_read_hires_parent' loads parent orography data from netCDF @@ -815,7 +818,7 @@ subroutine mn_orog_read_hires_parent(npx, npy, refine, pelist, surface_dir, filt real, allocatable, intent(out) :: orog_grid(:,:) !< Output orography grid real, allocatable, intent(out) :: orog_std_grid(:,:) !< Output orography standard deviation grid real, allocatable, intent(out) :: ls_mask_grid(:,:) !< Output land sea mask grid - real, allocatable, intent(out) :: land_frac_grid(:,:)!< Output land fraction grid + real(kind=kind_phys), allocatable, intent(out) :: land_frac_grid(:,:)!< Output land fraction grid integer, intent(in) :: parent_tile !< Parent tile number integer :: nx_cubic, nx, ny, fp_nx, fp_ny, mid_nx, mid_ny @@ -857,6 +860,110 @@ subroutine mn_orog_read_hires_parent(npx, npy, refine, pelist, surface_dir, filt end subroutine mn_orog_read_hires_parent + !>@brief The subroutine 'mn_replace_low_values' replaces low values with a default value. + subroutine mn_replace_low_values(data_grid, low_value, new_value) + real, _ALLOCATABLE, intent(inout) :: data_grid(:,:) !< 2D grid of data + real, intent(in) :: low_value !< Low value to check for; e.g. negative or fill value + real, intent(in) :: new_value !< Value to replace low value with + + integer :: i, j + + do i=lbound(data_grid,1),ubound(data_grid,1) + do j=lbound(data_grid,2),ubound(data_grid,2) + if (data_grid(i,j) .le. low_value) data_grid(i,j) = new_value + enddo + enddo + end subroutine mn_replace_low_values + + subroutine mn_static_read_ls(static_ls, npx, npy, refine, pelist, surface_dir, tile_num, terrain_smoother, filtered_terrain) + type(mn_land_mask_grids), intent(inout) :: static_ls + integer, intent(in) :: npx, npy, refine, tile_num !< Number of x,y points and nest refinement, (parent) tile number + integer, allocatable, intent(in) :: pelist(:) !< PE list for fms2_io + character(len=*), intent(in) :: surface_dir !< Surface directory + integer, intent(in) :: terrain_smoother + logical, intent(in) :: filtered_terrain + + ! If terrain_smoother method 1 is chosen, we need the parent coarse terrain + if (terrain_smoother .eq. 1) then + if (filtered_terrain) then + call mn_static_read_hires(npx, npy, refine, pelist, surface_dir, "oro_data", "orog_filt", static_ls%orog_grid, tile_num) + else + call mn_static_read_hires(npx, npy, refine, pelist, surface_dir, "oro_data", "orog_raw", static_ls%orog_grid, tile_num) + endif + endif + + ! Read in coarse resolution land sea mask to use for masked interpolations; factor in lakes as well + call mn_static_read_hires(npx, npy, refine, pelist, surface_dir, "oro_data", "slmsk", static_ls%ls_mask_grid, tile_num) + call mn_static_read_hires(npx, npy, refine, pelist, surface_dir, "oro_data", "land_frac", static_ls%land_frac_grid, tile_num) + + !! Lat lons for debugging + call mn_static_read_hires(npx, npy, refine, pelist, trim(surface_dir), "oro_data", "geolat", static_ls%geolat_grid, tile_num) + call mn_static_read_hires(npx, npy, refine, pelist, trim(surface_dir), "oro_data", "geolon", static_ls%geolon_grid, tile_num) + + ! Need parent soil type to determine lakes + call mn_static_read_hires(npx, npy, refine, pelist, trim(surface_dir), "soil_type", "soil_type", static_ls%soil_type_grid, tile_num) + ! To match initialization behavior, set any -999s to 0 in soil_type + call mn_replace_low_values(static_ls%soil_type_grid, -100.0, 0.0) + + end subroutine mn_static_read_ls + + subroutine mn_static_read_fix(static_fix, npx, npy, refine, pelist, surface_dir, tile_num, month) + type(mn_fix_grids), intent(inout) :: static_fix + integer, allocatable, intent(in) :: pelist(:) !< PE list for fms2_io + character(len=*), intent(in) :: surface_dir !< Surface directory + integer, intent(in) :: npx, npy, refine, tile_num, month !< Number of x,y points and nest refinement, (parent) tile number + + call mn_static_read_hires(npx, npy, refine, pelist, surface_dir, "substrate_temperature", "substrate_temperature", static_fix%deep_soil_temp_grid, tile_num) + ! set any -999s to +4C + call mn_replace_low_values(static_fix%deep_soil_temp_grid, -100.0, 277.0) + + + !! TODO investigate reading high-resolution veg_frac and veg_greenness + !call mn_static_read_hires(npx, npy, refine, trim(Moving_nest(child_grid_num)%mn_flag%surface_dir), "", mn_static%veg_frac_grid) + + call mn_static_read_hires(npx, npy, refine, pelist, trim(surface_dir), "vegetation_type", "vegetation_type", static_fix%veg_type_grid, tile_num) + ! To match initialization behavior, set any -999s to 0 in veg_type + call mn_replace_low_values(static_fix%veg_type_grid, -100.0, 0.0) + + + call mn_static_read_hires(npx, npy, refine, pelist, trim(surface_dir), "slope_type", "slope_type", static_fix%slope_type_grid, tile_num) + ! To match initialization behavior, set any -999s to 0 in slope_type + call mn_replace_low_values(static_fix%slope_type_grid, -100.0, 0.0) + + + call mn_static_read_hires(npx, npy, refine, pelist, trim(surface_dir), "maximum_snow_albedo", "maximum_snow_albedo", static_fix%max_snow_alb_grid, tile_num) + ! Set any -999s to 0.5 + call mn_replace_low_values(static_fix%max_snow_alb_grid, -100.0, 0.5) + + ! Albedo fraction -- read and calculate + call mn_static_read_hires(npx, npy, refine, pelist, trim(surface_dir), "facsf", "facsf", static_fix%facsf_grid, tile_num) + + call alloc_set_facwf(static_fix) + + ! Additional albedo variables + ! black sky = strong cosz -- direct sunlight + ! white sky = weak cosz -- diffuse light + + ! alvsf = visible strong cosz = visible_black_sky_albedo + ! alvwf = visible weak cosz = visible_white_sky_albedo + ! alnsf = near IR strong cosz = near_IR_black_sky_albedo + ! alnwf = near IR weak cosz = near_IR_white_sky_albedo + + call mn_static_read_hires(npx, npy, refine, pelist, trim(surface_dir), "snowfree_albedo", "visible_black_sky_albedo", static_fix%alvsf_grid, tile_num, time=month) + call mn_static_read_hires(npx, npy, refine, pelist, trim(surface_dir), "snowfree_albedo", "visible_white_sky_albedo", static_fix%alvwf_grid, tile_num, time=month) + + call mn_static_read_hires(npx, npy, refine, pelist, trim(surface_dir), "snowfree_albedo", "near_IR_black_sky_albedo", static_fix%alnsf_grid, tile_num, time=month) + call mn_static_read_hires(npx, npy, refine, pelist, trim(surface_dir), "snowfree_albedo", "near_IR_white_sky_albedo", static_fix%alnwf_grid, tile_num, time=month) + + ! Set the -999s to small value of 0.06, matching initialization code in chgres + + call mn_replace_low_values(static_fix%alvsf_grid, -100.0, 0.06) + call mn_replace_low_values(static_fix%alvwf_grid, -100.0, 0.06) + call mn_replace_low_values(static_fix%alnsf_grid, -100.0, 0.06) + call mn_replace_low_values(static_fix%alnwf_grid, -100.0, 0.06) + + end subroutine mn_static_read_fix + !>@brief The subroutine 'mn_static_read_hires_r4' loads high resolution data from netCDF !>@details Gathers a single variable from the netCDF file subroutine mn_static_read_hires_r4(npx, npy, refine, pelist, surface_dir, file_prefix, var_name, data_grid, parent_tile, time) @@ -1252,9 +1359,9 @@ subroutine mn_var_shift_data_r8_2d(data_var, interp_type, wt, ind, delta_i_c, de end subroutine mn_var_shift_data_r8_2d - !>@brief The subroutine 'mn_prog_shift_data_r4_3d' shifts the data for a variable on each nest PE + !>@brief The subroutine 'mn_prog_shift_data_r4_3d_highz' shifts the data for a variable on each nest PE !>@details For one single precision 3D variable - subroutine mn_var_shift_data_r4_3d(data_var, interp_type, wt, ind, delta_i_c, delta_j_c, x_refine, y_refine, is_fine_pe, nest_domain, position, nz) + subroutine mn_var_shift_data_r4_3d_highz(data_var, interp_type, wt, ind, delta_i_c, delta_j_c, x_refine, y_refine, is_fine_pe, nest_domain, position, nz) real*4, allocatable, intent(inout) :: data_var(:,:,:) !< Data variable integer, intent(in) :: interp_type !< Interpolation stagger type @@ -1265,6 +1372,22 @@ subroutine mn_var_shift_data_r4_3d(data_var, interp_type, wt, ind, delta_i_c, de type(nest_domain_type), intent(inout) :: nest_domain !< Nest domain structure integer, intent(in) :: position, nz !< Grid offset, number of vertical levels + call mn_var_shift_data_r4_3d_lowhighz(data_var, interp_type, wt, ind, delta_i_c, delta_j_c, x_refine, y_refine, is_fine_pe, nest_domain, position, 1, nz) + + end subroutine mn_var_shift_data_r4_3d_highz + !>@brief The subroutine 'mn_prog_shift_data_r4_3d_lowhighz' shifts the data for a variable on each nest PE + !>@details For one single precision 3D variable + subroutine mn_var_shift_data_r4_3d_lowhighz(data_var, interp_type, wt, ind, delta_i_c, delta_j_c, x_refine, y_refine, is_fine_pe, nest_domain, position, z_low, z_high) + + real*4, allocatable, intent(inout) :: data_var(:,:,:) !< Data variable + integer, intent(in) :: interp_type !< Interpolation stagger type + real, allocatable, intent(in) :: wt(:,:,:) !< Interpolation weight array + integer, allocatable, intent(in) :: ind(:,:,:) !< Fine to coarse index array + integer, intent(in) :: delta_i_c, delta_j_c, x_refine, y_refine !< delta i,j for nest move. Nest refinement. + logical, intent(in) :: is_fine_pe !< Is nest PE? + type(nest_domain_type), intent(inout) :: nest_domain !< Nest domain structure + integer, intent(in) :: position, z_low, z_high !< Grid offset, number of vertical levels + real*4, dimension(:,:,:), allocatable :: nbuffer, sbuffer, ebuffer, wbuffer type(bbox) :: north_fine, north_coarse ! step 4 type(bbox) :: south_fine, south_coarse @@ -1278,10 +1401,10 @@ subroutine mn_var_shift_data_r4_3d(data_var, interp_type, wt, ind, delta_i_c, de !! !!=========================================================== - call alloc_halo_buffer(nbuffer, north_fine, north_coarse, nest_domain, NORTH, position, nz) - call alloc_halo_buffer(sbuffer, south_fine, south_coarse, nest_domain, SOUTH, position, nz) - call alloc_halo_buffer(ebuffer, east_fine, east_coarse, nest_domain, EAST, position, nz) - call alloc_halo_buffer(wbuffer, west_fine, west_coarse, nest_domain, WEST, position, nz) + call alloc_halo_buffer(nbuffer, north_fine, north_coarse, nest_domain, NORTH, position, z_low, z_high) + call alloc_halo_buffer(sbuffer, south_fine, south_coarse, nest_domain, SOUTH, position, z_low, z_high) + call alloc_halo_buffer(ebuffer, east_fine, east_coarse, nest_domain, EAST, position, z_low, z_high) + call alloc_halo_buffer(wbuffer, west_fine, west_coarse, nest_domain, WEST, position, z_low, z_high) !==================================================== @@ -1310,10 +1433,10 @@ subroutine mn_var_shift_data_r4_3d(data_var, interp_type, wt, ind, delta_i_c, de !! !!=========================================================== - call fill_nest_from_buffer(interp_type, data_var, nbuffer, north_fine, north_coarse, nz, NORTH, x_refine, y_refine, wt, ind) - call fill_nest_from_buffer(interp_type, data_var, sbuffer, south_fine, south_coarse, nz, SOUTH, x_refine, y_refine, wt, ind) - call fill_nest_from_buffer(interp_type, data_var, ebuffer, east_fine, east_coarse, nz, EAST, x_refine, y_refine, wt, ind) - call fill_nest_from_buffer(interp_type, data_var, wbuffer, west_fine, west_coarse, nz, WEST, x_refine, y_refine, wt, ind) + call fill_nest_from_buffer(interp_type, data_var, nbuffer, north_fine, north_coarse, z_low, z_high, NORTH, x_refine, y_refine, wt, ind) + call fill_nest_from_buffer(interp_type, data_var, sbuffer, south_fine, south_coarse, z_low, z_high, SOUTH, x_refine, y_refine, wt, ind) + call fill_nest_from_buffer(interp_type, data_var, ebuffer, east_fine, east_coarse, z_low, z_high, EAST, x_refine, y_refine, wt, ind) + call fill_nest_from_buffer(interp_type, data_var, wbuffer, west_fine, west_coarse, z_low, z_high, WEST, x_refine, y_refine, wt, ind) endif deallocate(nbuffer) @@ -1321,12 +1444,12 @@ subroutine mn_var_shift_data_r4_3d(data_var, interp_type, wt, ind, delta_i_c, de deallocate(ebuffer) deallocate(wbuffer) - end subroutine mn_var_shift_data_r4_3d + end subroutine mn_var_shift_data_r4_3d_lowhighz - !>@brief The subroutine 'mn_prog_shift_data_r8_3d' shifts the data for a variable on each nest PE + !>@brief The subroutine 'mn_prog_shift_data_r8_3d_highz' shifts the data for a variable on each nest PE !>@details For one double precision 3D variable - subroutine mn_var_shift_data_r8_3d(data_var, interp_type, wt, ind, delta_i_c, delta_j_c, x_refine, y_refine, is_fine_pe, nest_domain, position, nz) + subroutine mn_var_shift_data_r8_3d_highz(data_var, interp_type, wt, ind, delta_i_c, delta_j_c, x_refine, y_refine, is_fine_pe, nest_domain, position, nz) real*8, allocatable, intent(inout) :: data_var(:,:,:) !< Data variable integer, intent(in) :: interp_type !< Interpolation stagger type @@ -1337,6 +1460,23 @@ subroutine mn_var_shift_data_r8_3d(data_var, interp_type, wt, ind, delta_i_c, de type(nest_domain_type), intent(inout) :: nest_domain !< Nest domain structure integer, intent(in) :: position, nz !< Grid offset, number vertical levels + call mn_var_shift_data_r8_3d_lowhighz(data_var, interp_type, wt, ind, delta_i_c, delta_j_c, x_refine, y_refine, is_fine_pe, nest_domain, position, 1, nz) + + end subroutine mn_var_shift_data_r8_3d_highz + + !>@brief The subroutine 'mn_prog_shift_data_r8_3d' shifts the data for a variable on each nest PE + !>@details For one double precision 3D variable + subroutine mn_var_shift_data_r8_3d_lowhighz(data_var, interp_type, wt, ind, delta_i_c, delta_j_c, x_refine, y_refine, is_fine_pe, nest_domain, position, z_low, z_high) + + real*8, allocatable, intent(inout) :: data_var(:,:,:) !< Data variable + integer, intent(in) :: interp_type !< Interpolation stagger type + real, allocatable, intent(in) :: wt(:,:,:) !< Interpolation weight array + integer, allocatable, intent(in) :: ind(:,:,:) !< Fine to coarse index array + integer, intent(in) :: delta_i_c, delta_j_c, x_refine, y_refine !< delta i,j for nest move. Nest refinement. + logical, intent(in) :: is_fine_pe !< Is nest PE? + type(nest_domain_type), intent(inout) :: nest_domain !< Nest domain structure + integer, intent(in) :: position, z_low, z_high !< Grid offset, number vertical levels + real*8, dimension(:,:,:), allocatable :: nbuffer, sbuffer, ebuffer, wbuffer type(bbox) :: north_fine, north_coarse ! step 4 type(bbox) :: south_fine, south_coarse @@ -1350,10 +1490,10 @@ subroutine mn_var_shift_data_r8_3d(data_var, interp_type, wt, ind, delta_i_c, de !! !!=========================================================== - call alloc_halo_buffer(nbuffer, north_fine, north_coarse, nest_domain, NORTH, position, nz) - call alloc_halo_buffer(sbuffer, south_fine, south_coarse, nest_domain, SOUTH, position, nz) - call alloc_halo_buffer(ebuffer, east_fine, east_coarse, nest_domain, EAST, position, nz) - call alloc_halo_buffer(wbuffer, west_fine, west_coarse, nest_domain, WEST, position, nz) + call alloc_halo_buffer(nbuffer, north_fine, north_coarse, nest_domain, NORTH, position, z_low, z_high) + call alloc_halo_buffer(sbuffer, south_fine, south_coarse, nest_domain, SOUTH, position, z_low, z_high) + call alloc_halo_buffer(ebuffer, east_fine, east_coarse, nest_domain, EAST, position, z_low, z_high) + call alloc_halo_buffer(wbuffer, west_fine, west_coarse, nest_domain, WEST, position, z_low, z_high) !==================================================== ! Passes data from coarse grid to fine grid's halo buffers; requires nest_domain to be intent(inout) @@ -1380,10 +1520,10 @@ subroutine mn_var_shift_data_r8_3d(data_var, interp_type, wt, ind, delta_i_c, de !! !!=========================================================== - call fill_nest_from_buffer(interp_type, data_var, nbuffer, north_fine, north_coarse, nz, NORTH, x_refine, y_refine, wt, ind) - call fill_nest_from_buffer(interp_type, data_var, sbuffer, south_fine, south_coarse, nz, SOUTH, x_refine, y_refine, wt, ind) - call fill_nest_from_buffer(interp_type, data_var, ebuffer, east_fine, east_coarse, nz, EAST, x_refine, y_refine, wt, ind) - call fill_nest_from_buffer(interp_type, data_var, wbuffer, west_fine, west_coarse, nz, WEST, x_refine, y_refine, wt, ind) + call fill_nest_from_buffer(interp_type, data_var, nbuffer, north_fine, north_coarse, z_low, z_high, NORTH, x_refine, y_refine, wt, ind) + call fill_nest_from_buffer(interp_type, data_var, sbuffer, south_fine, south_coarse, z_low, z_high, SOUTH, x_refine, y_refine, wt, ind) + call fill_nest_from_buffer(interp_type, data_var, ebuffer, east_fine, east_coarse, z_low, z_high, EAST, x_refine, y_refine, wt, ind) + call fill_nest_from_buffer(interp_type, data_var, wbuffer, west_fine, west_coarse, z_low, z_high, WEST, x_refine, y_refine, wt, ind) endif deallocate(nbuffer) @@ -1391,7 +1531,7 @@ subroutine mn_var_shift_data_r8_3d(data_var, interp_type, wt, ind, delta_i_c, de deallocate(ebuffer) deallocate(wbuffer) - end subroutine mn_var_shift_data_r8_3d + end subroutine mn_var_shift_data_r8_3d_lowhighz !>@brief The subroutine 'mn_prog_shift_data_r4_4d' shifts the data for a variable on each nest PE @@ -2121,14 +2261,14 @@ subroutine mn_var_dump_3d_to_netcdf( data_var, is_fine_pe, domain_coarse, domain call output_grid_to_nc("GH", isd_fine, ied_fine, jsd_fine, jed_fine, nz, data_var, prefix_fine, var_name, time_step, domain_fine, position) else - if (this_tile == 6) then + !if (this_tile == 6) then !call mpp_get_compute_domain(domain_coarse, isc_coarse, iec_coarse, jsc_coarse, jec_coarse, position=position) call mpp_get_data_domain(domain_coarse, isd_coarse, ied_coarse, jsd_coarse, jed_coarse, position=position) !call mpp_get_memory_domain(domain_coarse, ism_coarse, iem_coarse, jsm_coarse, jem_coarse, position=position) call output_grid_to_nc("GH", isd_coarse, ied_coarse, jsd_coarse, jed_coarse, nz, data_var, prefix_coarse, var_name, time_step, domain_coarse, position) - endif + !endif endif end subroutine mn_var_dump_3d_to_netcdf @@ -2172,14 +2312,14 @@ subroutine mn_var_dump_2d_to_netcdf( data_var, is_fine_pe, domain_coarse, domain call output_grid_to_nc("GH", isd_fine, ied_fine, jsd_fine, jed_fine, data_var, prefix_fine, var_name, time_step, domain_fine, position) else - if (this_tile == 6) then + !if (this_tile == 6) then !call mpp_get_compute_domain(domain_coarse, isc_coarse, iec_coarse, jsc_coarse, jec_coarse, position=position) call mpp_get_data_domain(domain_coarse, isd_coarse, ied_coarse, jsd_coarse, jed_coarse, position=position) !call mpp_get_memory_domain(domain_coarse, ism_coarse, iem_coarse, jsm_coarse, jem_coarse, position=position) call output_grid_to_nc("GH", isd_coarse, ied_coarse, jsd_coarse, jed_coarse, data_var, prefix_coarse, var_name, time_step, domain_coarse, position) - endif + !endif endif end subroutine mn_var_dump_2d_to_netcdf @@ -2439,7 +2579,7 @@ subroutine assign_p_grids(parent_geo, p_grid, position) do i = 1, ubound(p_grid,1) ! centered grid version p_grid(i, j, :) = 0.0 - + if (2*i .gt. ubound(parent_geo%lons,1)) then print '("[ERROR] WDR PG_CLONi npe=",I0," 2*i=",I0," ubound=",I0)', mpp_pe(), 2*i, ubound(parent_geo%lons,1) elseif (2*j .gt. ubound(parent_geo%lons,2)) then @@ -2450,7 +2590,7 @@ subroutine assign_p_grids(parent_geo, p_grid, position) endif enddo enddo - + do i = 1, ubound(p_grid,1) do j = 1, ubound(p_grid,2) @@ -2465,7 +2605,7 @@ subroutine assign_p_grids(parent_geo, p_grid, position) if (p_grid(i,j,1) .ne. 0.0) num_vals = num_vals + 1 enddo enddo - + if (num_zeros .gt. 0) print '("[INFO] WDR set p_grid npe=",I0," num_zeros=",I0," full_zeros=",I0," num_vals=",I0" nxp=",I0," nyp=",I0," parent_geo%lats(",I0,",",I0,")"," p_grid(",I0,",",I0,",2)")', mpp_pe(), num_zeros, full_zeros, num_vals, parent_geo%nxp, parent_geo%nyp, ubound(parent_geo%lats,1), ubound(parent_geo%lats,2), ubound(p_grid,1), ubound(p_grid,2) @@ -2480,15 +2620,15 @@ subroutine assign_p_grids(parent_geo, p_grid, position) elseif (2*j-1 .gt. ubound(parent_geo%lons,2)) then print '("[ERROR] WDR PG_ULONj npe=",I0," 2*j-1=",I0," ubound=",I0)', mpp_pe(), 2*j-1, ubound(parent_geo%lons,2) else - + ! This seems correct p_grid(i, j, 1) = parent_geo%lons(2*i, 2*j-1) p_grid(i, j, 2) = parent_geo%lats(2*i, 2*j-1) endif enddo enddo - - + + do i = 1, ubound(p_grid,1) do j = 1, ubound(p_grid,2) @@ -2516,18 +2656,18 @@ subroutine assign_p_grids(parent_geo, p_grid, position) print '("[ERROR] WDR PG_VLONi npe=",I0," 2*i-1=",I0," ubound=",I0)', mpp_pe(), 2*i-1, ubound(parent_geo%lons,1) elseif (2*j .gt. ubound(parent_geo%lons,2)) then print '("[ERROR] WDR PG_VLONj npe=",I0," 2*j=",I0," ubound=",I0)', mpp_pe(), 2*j, ubound(parent_geo%lons,2) - else + else ! This seems correct p_grid(i, j, 1) = parent_geo%lons(2*i-1, 2*j) p_grid(i, j, 2) = parent_geo%lats(2*i-1, 2*j) endif enddo enddo - + do i = 1, ubound(p_grid,1) do j = 1, ubound(p_grid,2) - + if (p_grid(i,j,1) .eq. 0.0) then num_zeros = num_zeros + 1 if (p_grid(i,j,2) .eq. 0.0) then @@ -2538,13 +2678,13 @@ subroutine assign_p_grids(parent_geo, p_grid, position) if (p_grid(i,j,1) .ne. 0.0) num_vals = num_vals + 1 enddo enddo - + if (num_zeros .gt. 0) print '("[INFO] WDR set p_grid_v npe=",I0," num_zeros=",I0," full_zeros=",I0," num_vals=",I0" nxp=",I0," nyp=",I0," parent_geo%lats(",I0,",",I0,")"," p_grid(",I0,",",I0,",2)")', mpp_pe(), num_zeros, full_zeros, num_vals, parent_geo%nxp, parent_geo%nyp, ubound(parent_geo%lats,1), ubound(parent_geo%lats,2), ubound(p_grid,1), ubound(p_grid,2) - - - + + + endif - + end subroutine assign_p_grids @@ -2596,9 +2736,9 @@ subroutine calc_inside(p_grid, ic, jc, n_grid1, n_grid2, istag, jstag, is_inside real(kind=R_GRID), intent(in) :: n_grid1, n_grid2 integer, intent(in) :: ic, jc, istag, jstag logical, intent(out) :: is_inside - logical, intent(in) :: verbose + logical, intent(in) :: verbose + - real(kind=R_GRID) :: max1, max2, min1, min2, eps max1 = max(p_grid(ic,jc,1), p_grid(ic,jc+1,1), p_grid(ic,jc+1,1), p_grid(ic+1,jc+1,1), p_grid(ic+1,jc+1,1), p_grid(ic+1,jc,1)) @@ -2606,21 +2746,21 @@ subroutine calc_inside(p_grid, ic, jc, n_grid1, n_grid2, istag, jstag, is_inside min1 = min(p_grid(ic,jc,1), p_grid(ic,jc+1,1), p_grid(ic,jc+1,1), p_grid(ic+1,jc+1,1), p_grid(ic+1,jc+1,1), p_grid(ic+1,jc,1)) min2 = min(p_grid(ic,jc,2), p_grid(ic,jc+1,2), p_grid(ic,jc+1,2), p_grid(ic+1,jc+1,2), p_grid(ic+1,jc+1,2), p_grid(ic+1,jc,2)) - + is_inside = .False. - + eps = 0.00001 !eps = 0.000001 if (n_grid1 .le. max1+eps .and. n_grid1 .ge. min1-eps) then if (n_grid2 .le. max2+eps .and. n_grid2 .ge. min2-eps) then is_inside = .True. - !if (verbose) print '("[INFO] WDR is_inside TRUE npe=",I0," ic=",I0," jc=",I0," n_grid1=",F12.8," min1=",F12.8," max1=",F12.8," n_grid2=",F12.8," min2=",F12.8," max2=", F12.8," p_grid(",I0,",",I0,",2) istag=",I0," jstag=",I0)', mpp_pe(), ic, jc, n_grid1, min1, max1, n_grid2, min2, max2, ubound(p_grid,1), ubound(p_grid,2), istag, jstag +! if (verbose) print '("[INFO] WDR is_inside TRUE npe=",I0," ic=",I0," jc=",I0," n_grid1=",F12.8," min1=",F12.8," max1=",F12.8," n_grid2=",F12.8," min2=",F12.8," max2=", F12.8," p_grid(",I0,",",I0,",2) istag=",I0," jstag=",I0)', mpp_pe(), ic, jc, n_grid1, min1, max1, n_grid2, min2, max2, ubound(p_grid,1), ubound(p_grid,2), istag, jstag endif endif if (verbose .and. .not. is_inside) then - print '("[INFO] WDR is_inside FALSE npe=",I0," ic=",I0," jc=",I0," n_grid1=",F12.8," min1=",F12.8," max1=",F12.8," n_grid2=",F12.8," min2=",F12.8," max2=", F10.4," p_grid(",I0,",",I0,",2) istag=",I0," jstag=",I0)', mpp_pe(), ic, jc, n_grid1, min1, max1, n_grid2, min2, max2, ubound(p_grid,1), ubound(p_grid,2), istag, jstag + print '("[WARN] WDR is_inside FALSE npe=",I0," ic=",I0," jc=",I0," n_grid1=",F12.8," min1=",F12.8," max1=",F12.8," n_grid2=",F12.8," min2=",F12.8," max2=", F10.4," p_grid(",I0,",",I0,",2) istag=",I0," jstag=",I0)', mpp_pe(), ic, jc, n_grid1, min1, max1, n_grid2, min2, max2, ubound(p_grid,1), ubound(p_grid,2), istag, jstag endif @@ -2672,7 +2812,7 @@ subroutine calc_nest_halo_weights(bbox_fine, bbox_coarse, p_grid, n_grid, wt, is ic = ind(i,j,1) if (ic+1 .gt. ubound(p_grid, 1)) print '("[ERROR] WDR CALCWT off end of p_grid i npe=",I0," ic+1=",I0," bound=",I0)', mpp_pe(), ic+1, ubound(p_grid,1) if (jc+1 .gt. ubound(p_grid, 2)) print '("[ERROR] WDR CALCWT off end of p_grid j npe=",I0," jc+1=",I0," bound=",I0)', mpp_pe(), jc+1, ubound(p_grid,2) - + ! dist2side_latlon takes points in longitude-latitude coordinates. dist1 = dist2side_latlon(p_grid(ic,jc,:), p_grid(ic,jc+1,:), n_grid(i,j,:)) dist2 = dist2side_latlon(p_grid(ic,jc+1,:), p_grid(ic+1,jc+1,:), n_grid(i,j,:)) @@ -2681,33 +2821,33 @@ subroutine calc_nest_halo_weights(bbox_fine, bbox_coarse, p_grid, n_grid, wt, is call calc_inside(p_grid, ic, jc, n_grid(i,j,1), n_grid(i,j,2), istag, jstag, is_inside, .True.) -! if (.not. is_inside) then -! adjusted = .False. -! -! do di = -2,2 -! do dj = -2,1 -! if (.not. adjusted) then -! call calc_inside(p_grid, ic+di, jc+dj, n_grid(i,j,1), n_grid(i,j,2), istag, jstag, is_inside, .False.) -! if (is_inside) then -! ic = ic + di -! jc = jc + dj -! -! print '("[INFO] WDR is_inside UPDATED npe=",I0," ic=",I0," jc=",I0," istart_coarse=",I0," jstart_coarse=",I0," i=",I0," j=",I0," di=",I0," dj=",I0," n_grid1=",F12.8," n_grid2=",F12.8," istag=",I0," jstag=",I0)', mpp_pe(), ic, jc, istart_coarse, jstart_coarse, i, j, di, dj, n_grid(i,j,1), n_grid(i,j,2), istag, jstag - -! dist1 = dist2side_latlon(p_grid(ic,jc,:), p_grid(ic,jc+1,:), n_grid(i,j,:)) -! dist2 = dist2side_latlon(p_grid(ic,jc+1,:), p_grid(ic+1,jc+1,:), n_grid(i,j,:)) -! dist3 = dist2side_latlon(p_grid(ic+1,jc+1,:), p_grid(ic+1,jc,:), n_grid(i,j,:)) -! dist4 = dist2side_latlon(p_grid(ic,jc,:), p_grid(ic+1,jc,:), n_grid(i,j,:)) -! -! adjusted = .True. -! endif -! endif -! enddo -! enddo -! if (.not. adjusted) print '("[ERROR] WDR is_inside UPDATE FAILED npe=",I0," i=",I0," j=",I0," ic=",I0," jc=",I0," n_grid1=",F12.8," n_grid2=",F12.8," istag=",I0," jstag=",I0)', mpp_pe(), i, j, ic, jc, n_grid(i,j,1), n_grid(i,j,2), istag, jstag -! -! endif - + if (.not. is_inside) then + adjusted = .False. + + do di = -2,2 + do dj = -2,1 + if (.not. adjusted) then + call calc_inside(p_grid, ic+di, jc+dj, n_grid(i,j,1), n_grid(i,j,2), istag, jstag, is_inside, .False.) + if (is_inside) then + ic = ic + di + jc = jc + dj + + print '("[INFO] WDR is_inside UPDATED npe=",I0," ic=",I0," jc=",I0," istart_coarse=",I0," jstart_coarse=",I0," i=",I0," j=",I0," di=",I0," dj=",I0," n_grid1=",F12.8," n_grid2=",F12.8," istag=",I0," jstag=",I0)', mpp_pe(), ic, jc, istart_coarse, jstart_coarse, i, j, di, dj, n_grid(i,j,1), n_grid(i,j,2), istag, jstag + + dist1 = dist2side_latlon(p_grid(ic,jc,:), p_grid(ic,jc+1,:), n_grid(i,j,:)) + dist2 = dist2side_latlon(p_grid(ic,jc+1,:), p_grid(ic+1,jc+1,:), n_grid(i,j,:)) + dist3 = dist2side_latlon(p_grid(ic+1,jc+1,:), p_grid(ic+1,jc,:), n_grid(i,j,:)) + dist4 = dist2side_latlon(p_grid(ic,jc,:), p_grid(ic+1,jc,:), n_grid(i,j,:)) + + adjusted = .True. + endif + endif + enddo + enddo + if (.not. adjusted) print '("[ERROR] WDR is_inside UPDATE FAILED npe=",I0," i=",I0," j=",I0," ic=",I0," jc=",I0," n_grid1=",F12.8," n_grid2=",F12.8," istag=",I0," jstag=",I0)', mpp_pe(), i, j, ic, jc, n_grid(i,j,1), n_grid(i,j,2), istag, jstag + + endif + old_weight = wt(i,j,:) wt(i,j,1)=dist2*dist3 ! ic, jc weight @@ -2723,9 +2863,9 @@ subroutine calc_nest_halo_weights(bbox_fine, bbox_coarse, p_grid, n_grid, wt, is call mpp_error(WARNING, "WARNING: calc_nest_halo_weights sum of weights is zero.") wt(i,j,:) = 0.25 - + else - wt(i,j,:)=wt(i,j,:)/sum + wt(i,j,:)=wt(i,j,:)/sum endif diff_weight = old_weight - wt(i,j,:) @@ -2733,7 +2873,7 @@ subroutine calc_nest_halo_weights(bbox_fine, bbox_coarse, p_grid, n_grid, wt, is if (abs(diff_weight(k)) .ge. 0.01) then print '("[WARN] WDR DIFFWT npe=",I0," old_wt=",F10.6," wt(",I0,",",I0,",",I0,")=",F10.6," diff=",F10.6," istag=",I0," jstag=",I0)', & mpp_pe(), old_weight(k), i, j, k, wt(i,j,k), diff_weight(k), istag, jstag - + endif enddo enddo @@ -2743,4 +2883,3 @@ subroutine calc_nest_halo_weights(bbox_fine, bbox_coarse, p_grid, n_grid, wt, is end subroutine calc_nest_halo_weights end module fv_moving_nest_mod - diff --git a/moving_nest/fv_moving_nest_main.F90 b/moving_nest/fv_moving_nest_main.F90 index 34af608c2..c4f53bc2c 100644 --- a/moving_nest/fv_moving_nest_main.F90 +++ b/moving_nest/fv_moving_nest_main.F90 @@ -106,13 +106,17 @@ module fv_moving_nest_main_mod use fv_moving_nest_types_mod, only: allocate_fv_moving_nest_prog_type, allocate_fv_moving_nest_physics_type use fv_moving_nest_types_mod, only: deallocate_fv_moving_nests use fv_moving_nest_types_mod, only: Moving_nest + use fv_moving_nest_types_mod, only: mn_apply_lakes, mn_overwrite_with_nest_init_values, alloc_set_facwf + use fv_moving_nest_types_mod, only: mn_static_overwrite_ls_from_nest, mn_static_overwrite_fix_from_nest + use fv_moving_nest_types_mod, only: deallocate_land_mask_grids, deallocate_fix_grids ! Prognostic variable routines use fv_moving_nest_mod, only: mn_prog_fill_intern_nest_halos, mn_prog_fill_nest_halos_from_parent, & mn_prog_dump_to_netcdf, mn_prog_shift_data + use fv_moving_nest_mod, only: mn_static_read_ls, mn_static_read_fix ! Physics variable routines use fv_moving_nest_physics_mod, only: mn_phys_fill_intern_nest_halos, mn_phys_fill_nest_halos_from_parent, & - mn_phys_dump_to_netcdf, mn_phys_shift_data, mn_phys_reset_sfc_props, move_nsst + mn_phys_dump_to_netcdf, mn_phys_shift_data, mn_phys_reset_sfc_props, move_nsst, mn_phys_set_slmsk ! Metadata routines use fv_moving_nest_mod, only: mn_meta_move_nest, mn_meta_recalc, mn_meta_reset_gridstruct, mn_shift_index @@ -165,6 +169,9 @@ module fv_moving_nest_main_mod integer :: id_movnestTot integer, save :: output_step = 0 + type(mn_surface_grids), save :: mn_static + + contains !>@brief The subroutine 'update_moving_nest' decides whether the nest should be moved, and if so, performs the move. @@ -247,6 +254,236 @@ end subroutine nest_tracker_end + + subroutine log_landsea_mask(Atm_block, IPD_control, IPD_data, time_step, parent_grid_num, child_grid_num) + type(block_control_type), intent(in) :: Atm_block !< Physics block layout + type(IPD_control_type), intent(in) :: IPD_control !< Physics metadata + type(IPD_data_type), intent(in) :: IPD_data(:) !< Physics variable data + type(time_type), intent(in) :: time_step !< Current timestep + integer, intent(in) :: parent_grid_num, child_grid_num + + + character(len=160) :: line + character(len=1) :: mask_char + character(len=1) :: num_char + integer :: i,j + integer :: nb, blen, ix, i_pe, j_pe, i_idx, j_idx, refine + integer :: ioffset, joffset + real :: local_slmsk(Atm(2)%bd%isd:Atm(2)%bd%ied, Atm(2)%bd%jsd:Atm(2)%bd%jed) + integer :: nz, this_pe, n + integer :: num_land, num_water + + this_pe = mpp_pe() + n = mygrid + + refine = Atm(child_grid_num)%neststruct%refinement + ioffset = Atm(child_grid_num)%neststruct%ioffset + joffset = Atm(child_grid_num)%neststruct%joffset + + do i=lbound(Atm(n)%oro,1), ubound(Atm(n)%oro,1) + line = "" + do j=lbound(Atm(n)%oro,2), ubound(Atm(n)%oro,2) + !print '("[INFO] WDR oro size npe=",I0," is_allocated=",L1)', this_pe, allocated(Atm(n)%oro) + !print '("[INFO] WDR oro size npe=",I0," i=",I0,"-",I0," j=",I0,"-",I0)', this_pe, lbound(Atm(n)%oro,1), ubound(Atm(n)%oro,1), lbound(Atm(n)%oro,2), ubound(Atm(n)%oro,2) + if (Atm(n)%oro(i,j) .eq. 1) then + ! land + line = trim(line) // "+" + elseif (Atm(n)%oro(i,j) .eq. 2) then + ! Water + line = trim(line) // "." + else + ! Unknown + line = trim(line) // "X" + endif + enddo + !print '("[INFO] WDR oro npe=",I0," time=",I0," i=",I0," ",A80)',this_pe,a_step,i,trim(line) + + enddo + + + local_slmsk = 8 + !print '("[INFO] WDR local_slmsk size npe=",I0," i=",I0,"-",I0," j=",I0,"-",I0," n=",I0)', this_pe, lbound(local_slmsk,1), ubound(local_slmsk,1), lbound(local_slmsk,2), ubound(local_slmsk,2), n + + do nb = 1,Atm_block%nblks + blen = Atm_block%blksz(nb) + do ix = 1, blen + i_pe = Atm_block%index(nb)%ii(ix) + j_pe = Atm_block%index(nb)%jj(ix) + + !print '("[INFO] WDR local_slmsk npe=",I0," i_pe=",I0," j_pe=",I0)', this_pe, i_pe, j_pe + + local_slmsk(i_pe, j_pe) = IPD_data(nb)%Sfcprop%slmsk(ix) + + if (allocated(Moving_nest)) then + if (allocated(Moving_nest(n)%mn_phys%slmsk)) then + if (int(local_slmsk(i_pe,j_pe)) .ne. 8) then + if (int(local_slmsk(i_pe,j_pe)) .ne. int(Moving_nest(n)%mn_phys%slmsk(i_pe,j_pe))) then + print '("[INFO] WDR mismatch local_slmsk_lake npe=",I0," time=",I3," i_pe=",I3," j_pe=",I3," slmsk=",I0," phys%slmsk=",I0," soil_type_grid=",I0," phys%soil_type=",I0," ipd%landfrac=",F10.5," land_frac_grid=",F12.5," ipd%lakefrac=",F10.5," ipd%oceanfrac=",F10.5)', & + this_pe,a_step,i_pe,j_pe, int(local_slmsk(i_pe,j_pe)), & + int(Moving_nest(n)%mn_phys%slmsk(i_pe,j_pe)), & + int(IPD_data(nb)%Sfcprop%stype(ix)), & + int(mn_static%fp_ls%soil_type_grid((ioffset-1)*refine+i_pe, (joffset-1)*refine+j_pe)), & + IPD_data(nb)%Sfcprop%landfrac(ix), & + int(mn_static%fp_ls%land_frac_grid((ioffset-1)*refine+i_pe, (joffset-1)*refine+j_pe)), & + IPD_data(nb)%Sfcprop%lakefrac(ix), & + IPD_data(nb)%Sfcprop%oceanfrac(ix) + endif + endif + endif + endif + enddo + enddo + + print '("[INFO] WDR local_slmsk size npe=",I0," i=",I0,"-",I0," j=",I0,"-",I0)', this_pe, lbound(local_slmsk,1), ubound(local_slmsk,1), lbound(local_slmsk,2), ubound(local_slmsk,2) + + line = "" + do j=lbound(local_slmsk,2), ubound(local_slmsk,2) + write(num_char, "(I1)"), mod(j,10) + line = trim(line) // trim(num_char) + enddo + print '("[INFO] WDR local_slmsk_lake npe=",I0," time=",I3," i=",I3," ",A60)',this_pe,a_step,-99,trim(line) + + do i=lbound(local_slmsk,1), ubound(local_slmsk,1) + line = "" + num_land = 0 + num_water = 0 + + do j=lbound(local_slmsk,2), ubound(local_slmsk,2) + + if (local_slmsk(i,j) .eq. 1) then + ! land + line = trim(line) // "+" + num_land = num_land + 1 + elseif (local_slmsk(i,j) .eq. 2) then + ! Water + line = trim(line) // "T" + elseif (local_slmsk(i,j) .eq. 0) then + ! Zero == lake? + line = trim(line) // "." + num_water = num_water + 1 + elseif (local_slmsk(i,j) .eq. 8) then + ! Missing/edge + line = trim(line) // "M" + else + ! Unknown + print '("[INFO] WDR local_slmsk_lake npe=",I0," time=",I3," i=",I3," j=",I3," slmsk=",E12.5)',this_pe,a_step,i,j, local_slmsk(i,j) + write (mask_char, "(I1)") int(local_slmsk(i,j)) + line = trim(line) // mask_char + endif + enddo + print '("[INFO] WDR local_slmsk_lake npe=",I0," time=",I3," i=",I3," ",A60," ",I2," ",I2)',this_pe,a_step,i,trim(line), num_land, num_water + enddo + end subroutine log_landsea_mask + + + subroutine validate_geo_coords(tag, geo_grid, nest_geo_grid, refine, ioffset, joffset) + character(len=*) :: tag + real(kind=kind_phys), allocatable, intent(in) :: geo_grid(:,:) + real(kind=kind_phys), allocatable, intent(in) :: nest_geo_grid(:,:) + integer, intent(in) :: refine, ioffset, joffset + + integer :: i,j, this_pe + real(kind=kind_phys) :: diff + + this_pe = mpp_pe() + + do i = lbound(nest_geo_grid,1), ubound(nest_geo_grid,1) + do j = lbound(nest_geo_grid,2), ubound(nest_geo_grid,2) + + diff = nest_geo_grid(i,j) - geo_grid((ioffset-1)*refine+i, (joffset-1)*refine+j) + print '("[INFO] WDR VALIDATEGEO tag=",A3," npe=",I0," i=",I0," j=",I0," diff=",F20.12," nest_val=",F16.12)', tag, this_pe, i, j, diff, nest_geo_grid(i,j) + + enddo + enddo + + end subroutine validate_geo_coords + + + + subroutine validate_navigation_fields(tag, Atm_block, IPD_control, IPD_data, parent_grid_num, child_grid_num) + character(len=*) :: tag + type(block_control_type), intent(in) :: Atm_block !< Physics block layout + type(IPD_control_type), intent(in) :: IPD_control !< Physics metadata + type(IPD_data_type), intent(in) :: IPD_data(:) !< Physics variable data + integer, intent(in) :: parent_grid_num, child_grid_num + + + character(len=160) :: line + character(len=1) :: mask_char + character(len=1) :: num_char + integer :: i,j + integer :: nb, blen, ix, i_pe, j_pe, i_idx, j_idx, refine + integer :: ioffset, joffset + real :: local_slmsk(Atm(2)%bd%isd:Atm(2)%bd%ied, Atm(2)%bd%jsd:Atm(2)%bd%jed) + integer :: nz, this_pe, n + integer :: num_land, num_water + + this_pe = mpp_pe() + n = mygrid + + refine = Atm(child_grid_num)%neststruct%refinement + ioffset = Atm(child_grid_num)%neststruct%ioffset + joffset = Atm(child_grid_num)%neststruct%joffset + + local_slmsk = 8 + print '("[INFO] WDR VALIDATE local_slmsk size npe=",I0," i=",I0,"-",I0," j=",I0,"-",I0," n=",I0)', this_pe, lbound(local_slmsk,1), ubound(local_slmsk,1), lbound(local_slmsk,2), ubound(local_slmsk,2), n + + do nb = 1,Atm_block%nblks + blen = Atm_block%blksz(nb) + do ix = 1, blen + i_pe = Atm_block%index(nb)%ii(ix) + j_pe = Atm_block%index(nb)%jj(ix) + + !print '("[INFO] WDR local_slmsk npe=",I0," i_pe=",I0," j_pe=",I0)', this_pe, i_pe, j_pe + + local_slmsk(i_pe, j_pe) = IPD_data(nb)%Sfcprop%slmsk(ix) + + if (allocated(Moving_nest)) then + if (allocated(Moving_nest(n)%mn_phys%slmsk)) then + !print '("[INFO] WDR VALIDATE local_slmsk npe=",I0," i_pe=",I0," j_pe=",I0)', this_pe, i_pe, j_pe + + if (int(local_slmsk(i_pe,j_pe)) .ne. 8) then + if (int(local_slmsk(i_pe,j_pe)) .ne. int(Moving_nest(n)%mn_phys%slmsk(i_pe,j_pe))) then + print '("[INFO] WDR mismatch VALIDATE A tag=",A4," npe=",I0," time=",I3," i_pe=",I3," j_pe=",I3," IPD%slmsk=",I0," phys%slmsk=",I0," fp_slmsk=",I0," soil_type_grid=",I0," phys%soil_type=",I0," ipd%landfrac=",F10.5," land_frac_grid=",F12.5," ipd%lakefrac=",F10.5," ipd%oceanfrac=",F10.5)', & + tag, this_pe,a_step,i_pe,j_pe, int(local_slmsk(i_pe,j_pe)), & + int(Moving_nest(n)%mn_phys%slmsk(i_pe,j_pe)), & + int(mn_static%fp_ls%ls_mask_grid((ioffset-1)*refine+i_pe, (joffset-1)*refine+j_pe)), & + int(IPD_data(nb)%Sfcprop%stype(ix)), & + int(mn_static%fp_ls%soil_type_grid((ioffset-1)*refine+i_pe, (joffset-1)*refine+j_pe)), & + IPD_data(nb)%Sfcprop%landfrac(ix), & + int(mn_static%fp_ls%land_frac_grid((ioffset-1)*refine+i_pe, (joffset-1)*refine+j_pe)), & + IPD_data(nb)%Sfcprop%lakefrac(ix), & + IPD_data(nb)%Sfcprop%oceanfrac(ix) + endif + + +! if ((i_pe .eq. 149 .and. j_pe .eq. 169) .or.(i_pe .eq. 152 .and. j_pe .eq. 169) .or. int(local_slmsk(i_pe,j_pe)) .ne. int(mn_static%ls_mask_grid((ioffset-1)*refine+i_pe, (joffset-1)*refine+j_pe))) then + if (int(local_slmsk(i_pe,j_pe)) .ne. int(mn_static%fp_ls%ls_mask_grid((ioffset-1)*refine+i_pe, (joffset-1)*refine+j_pe))) then + print '("[INFO] WDR mismatch VALIDATE B tag=",A4," npe=",I0," time=",I3," i_pe=",I3," j_pe=",I3," IPD%slmsk=",I0," phys%slmsk=",I0," fp_slmsk=",I0," soil_type_grid=",I0," phys%soil_type=",I0," ipd%landfrac=",F10.5," land_frac_grid=",F12.5," ipd%lakefrac=",F10.5," ipd%oceanfrac=",F10.5)', & + tag, this_pe,a_step,i_pe,j_pe, int(local_slmsk(i_pe,j_pe)), & + int(Moving_nest(n)%mn_phys%slmsk(i_pe,j_pe)), & + int(mn_static%fp_ls%ls_mask_grid((ioffset-1)*refine+i_pe, (joffset-1)*refine+j_pe)), & + int(IPD_data(nb)%Sfcprop%stype(ix)), & + int(mn_static%fp_ls%soil_type_grid((ioffset-1)*refine+i_pe, (joffset-1)*refine+j_pe)), & + IPD_data(nb)%Sfcprop%landfrac(ix), & + int(mn_static%fp_ls%land_frac_grid((ioffset-1)*refine+i_pe, (joffset-1)*refine+j_pe)), & + IPD_data(nb)%Sfcprop%lakefrac(ix), & + IPD_data(nb)%Sfcprop%oceanfrac(ix) + endif + + + + endif + endif + endif + enddo + enddo + + + + end subroutine validate_navigation_fields + + !>@brief The subroutine 'dump_moving_nest' outputs native grid format data to netCDF files !>@details This subroutine exports model variables using FMS IO to netCDF files if tsvar_out is set to .True. subroutine dump_moving_nest(Atm_block, IPD_control, IPD_data, time_step) @@ -275,6 +512,18 @@ subroutine dump_moving_nest(Atm_block, IPD_control, IPD_data, time_step) ! if (tsvar_out) call mn_prog_dump_to_netcdf(Atm(n), a_step, "tsavar", is_fine_pe, domain_coarse, domain_fine, nz) ! if (tsvar_out) call mn_phys_dump_to_netcdf(Atm(n), Atm_block, IPD_control, IPD_data, a_step, "tsavar", is_fine_pe, domain_coarse, domain_fine, nz) !endif + !if (a_step .ge. 310) then + !if (mod(a_step, 80) .eq. 0 ) then + ! if (tsvar_out) call mn_phys_dump_to_netcdf(Atm(n), Atm_block, IPD_control, IPD_data, a_step, "tsavar", is_fine_pe, domain_coarse, domain_fine, nz) + !endif + + ! if (is_fine_pe) then + ! call validate_navigation_fields("DUMP", Atm_block, IPD_control, IPD_data, parent_grid_num, child_grid_num) + ! endif + + !if (this_pe .eq. 88 .or. this_pe .eq. 89) then + ! call log_landsea_mask(Atm_block, IPD_control, IPD_data, time_step, parent_grid_num, child_grid_num) + !endif end subroutine dump_moving_nest @@ -499,7 +748,6 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, integer, pointer :: ioffset, joffset real, pointer, dimension(:,:,:) :: grid, agrid type(domain2d), pointer :: domain_coarse, domain_fine - real(kind=R_GRID), pointer, dimension(:,:,:,:) :: grid_global ! Constants for mpp calls integer :: position = CENTER @@ -509,23 +757,17 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, integer :: x_refine, y_refine ! Currently equal, but allows for future flexibility logical :: is_fine_pe - ! TODO read halo size from the namelist instead to allow nest refinement > 3 - integer :: ehalo = 3 - integer :: whalo = 3 - integer :: nhalo = 3 - integer :: shalo = 3 integer :: extra_halo = 0 ! Extra halo for moving nest routines integer :: istart_fine, iend_fine, jstart_fine, jend_fine integer :: istart_coarse, iend_coarse, jstart_coarse, jend_coarse integer :: nx, ny, nz, nx_cubic, ny_cubic - integer :: p_istart_fine, p_iend_fine, p_jstart_fine, p_jend_fine ! Parent tile data, saved between timesteps logical, save :: first_nest_move = .true. type(grid_geometry), save :: parent_geo type(grid_geometry), save :: fp_super_tile_geo - type(mn_surface_grids), save :: mn_static +! 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(:,:,:) @@ -541,9 +783,8 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, !real :: va(isd:ied,jsd:jed) logical :: filtered_terrain = .True. ! TODO set this from namelist - integer :: i, j, x, y, z, p, nn, n_moist + integer :: i, j integer :: parent_tile - logical :: found_nest_domain = .false. ! Variables to enable debugging use of mpp_sync logical :: debug_sync = .false. @@ -551,15 +792,12 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, integer :: pp, p1, p2 ! Variables for parent side of setup_aligned_nest() - integer :: isg, ieg, jsg, jeg, gid - integer :: isc_p, iec_p, jsc_p, jec_p - integer :: upoff, jind - integer :: ng, refinement integer :: npx, npy, npz, ncnst, pnats integer :: isc, iec, jsc, jec integer :: isd, ied, jsd, jed integer :: nq ! number of transported tracers integer :: is, ie, js, je, k ! For recalculation of omga + integer :: i_idx, j_idx integer, save :: output_step = 0 integer, allocatable :: pelist(:) character(len=16) :: errstring @@ -567,13 +805,21 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, integer :: year, month, day, hour, minute, second real(kind=R_GRID) :: pi = 4 * atan(1.0d0) real :: rad2deg + integer :: static_nest_num logical :: use_timers + real(kind=kind_phys):: maxSkinTempK rad2deg = 180.0 / pi - gid = mpp_pe() this_pe = mpp_pe() + ! Highest satellite observed skin temperatures on Earth are on the order of +70C/343K/+160F + ! Mildrexler, D. J., M. Zhao, and S. W. Running, 2011: Satellite Finds Highest Land Skin Temperatures on Earth. Bull. Amer. Meteor. Soc., 92, 855–860, + ! https://doi.org/10.1175/2011BAMS3067.1. + ! https://journals.ametsoc.org/view/journals/bams/92/7/2011bams3067_1.xml?tab_body=pdf + + maxSkinTempK = 273.15 + 80.0 + use_timers = Atm(n)%flagstruct%fv_timers allocate(pelist(mpp_npes())) @@ -608,9 +854,8 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, is_fine_pe = Atm(n)%neststruct%nested .and. ANY(Atm(n)%pelist(:) == this_pe) - if (first_nest_move) then - + call fv_moving_nest_init_clocks(Atm(n)%flagstruct%fv_timers) ! If NSST is turned off, do not move the NSST variables. @@ -621,12 +866,14 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, move_nsst=.true. endif + ! This will only allocate the mn_prog and mn_phys for the active Atm(n), not all of them ! The others can safely remain unallocated. call allocate_fv_moving_nest_prog_type(isd, ied, jsd, jed, npz, Moving_nest(n)%mn_prog) call allocate_fv_moving_nest_physics_type(isd, ied, jsd, jed, npz, move_physics, move_nsst, & - IPD_Control%lsoil, IPD_Control%nmtvr, IPD_Control%levs, IPD_Control%ntot2d, IPD_Control%ntot3d, & + IPD_Control%lsnow_lsm_lbound, IPD_Control%lsnow_lsm_ubound, IPD_Control%lsoil, & + IPD_Control%nmtvr, IPD_Control%levs, IPD_Control%ntot2d, IPD_Control%ntot3d, & Moving_nest(n)%mn_phys) endif @@ -740,90 +987,54 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, ! Also read in other static variables from the orography and surface files if (first_nest_move) then + ! TODO Compute this more flexibly for multiple moving nests + if (parent_tile .eq. 1) then + static_nest_num = 8 ! Regional + else + static_nest_num = 7 ! Global + endif + + !print '("[INFO] WDR NEST_NUM npe=",I0," is_regional=",L1," static_nest_num=",I0," parent_tile=",I0,", ntiles=",I0)', this_pe, Atm(n)%flagstruct%regional, static_nest_num, parent_tile, Atm(1)%flagstruct%ntiles ! TODO set pelist for the correct nest instead of hard-coded Atm(2)%pelist to allow multiple moving nests call mn_latlon_read_hires_parent(Atm(1)%npx, Atm(1)%npy, x_refine, Atm(2)%pelist, fp_super_tile_geo, & Moving_nest(child_grid_num)%mn_flag%surface_dir, parent_tile) - call mn_orog_read_hires_parent(Atm(1)%npx, Atm(1)%npy, x_refine, Atm(2)%pelist, & - Moving_nest(child_grid_num)%mn_flag%surface_dir, filtered_terrain, & - mn_static%orog_grid, mn_static%orog_std_grid, mn_static%ls_mask_grid, mn_static%land_frac_grid, parent_tile) - - ! If terrain_smoother method 1 is chosen, we need the parent coarse terrain - if (Moving_nest(n)%mn_flag%terrain_smoother .eq. 1) then - if (filtered_terrain) then - call mn_static_read_hires(Atm(1)%npx, Atm(1)%npy, 1, Atm(2)%pelist, Moving_nest(child_grid_num)%mn_flag%surface_dir, "oro_data", "orog_filt", mn_static%parent_orog_grid, parent_tile) - else - call mn_static_read_hires(Atm(1)%npx, Atm(1)%npy, 1, Atm(2)%pelist, Moving_nest(child_grid_num)%mn_flag%surface_dir, "oro_data", "orog_raw", mn_static%parent_orog_grid, parent_tile) - endif - endif - - call mn_static_read_hires(Atm(1)%npx, Atm(1)%npy, x_refine, Atm(2)%pelist, trim(Moving_nest(child_grid_num)%mn_flag%surface_dir), "substrate_temperature", "substrate_temperature", mn_static%deep_soil_temp_grid, parent_tile) - ! set any -999s to +4C - call mn_replace_low_values(mn_static%deep_soil_temp_grid, -100.0, 277.0) - - call mn_static_read_hires(Atm(1)%npx, Atm(1)%npy, x_refine, Atm(2)%pelist, trim(Moving_nest(child_grid_num)%mn_flag%surface_dir), "soil_type", "soil_type", mn_static%soil_type_grid, parent_tile) - ! To match initialization behavior, set any -999s to 0 in soil_type - call mn_replace_low_values(mn_static%soil_type_grid, -100.0, 0.0) - - - !! TODO investigate reading high-resolution veg_frac and veg_greenness - !call mn_static_read_hires(Atm(1)%npx, Atm(1)%npy, x_refine, trim(Moving_nest(child_grid_num)%mn_flag%surface_dir), "", mn_static%veg_frac_grid) + ! Read static parent land sea mask fields + call mn_static_read_ls(mn_static%parent_ls, Atm(1)%npx, Atm(1)%npy, 1, Atm(2)%pelist, trim(Moving_nest(child_grid_num)%mn_flag%surface_dir), parent_tile, Moving_nest(n)%mn_flag%terrain_smoother, filtered_terrain) - call mn_static_read_hires(Atm(1)%npx, Atm(1)%npy, x_refine, Atm(2)%pelist, trim(Moving_nest(child_grid_num)%mn_flag%surface_dir), "vegetation_type", "vegetation_type", mn_static%veg_type_grid, parent_tile) - ! To match initialization behavior, set any -999s to 0 in veg_type - call mn_replace_low_values(mn_static%veg_type_grid, -100.0, 0.0) + ! Read full panel + call mn_static_read_ls(mn_static%fp_ls, Atm(1)%npx, Atm(1)%npy, x_refine, Atm(2)%pelist, trim(Moving_nest(child_grid_num)%mn_flag%surface_dir), parent_tile, Moving_nest(n)%mn_flag%terrain_smoother, filtered_terrain) + ! Read static nest land sea mask fields + call mn_static_read_ls(mn_static%nest_ls, Atm(2)%npx, Atm(2)%npy, 1, Atm(2)%pelist, trim(Moving_nest(child_grid_num)%mn_flag%surface_dir) // "/..", static_nest_num, Moving_nest(n)%mn_flag%terrain_smoother, filtered_terrain) - call mn_static_read_hires(Atm(1)%npx, Atm(1)%npy, x_refine, Atm(2)%pelist, trim(Moving_nest(child_grid_num)%mn_flag%surface_dir), "slope_type", "slope_type", mn_static%slope_type_grid, parent_tile) - ! To match initialization behavior, set any -999s to 0 in slope_type - call mn_replace_low_values(mn_static%slope_type_grid, -100.0, 0.0) - - - call mn_static_read_hires(Atm(1)%npx, Atm(1)%npy, x_refine, Atm(2)%pelist, trim(Moving_nest(child_grid_num)%mn_flag%surface_dir), "maximum_snow_albedo", "maximum_snow_albedo", mn_static%max_snow_alb_grid, parent_tile) - ! Set any -999s to 0.5 - call mn_replace_low_values(mn_static%max_snow_alb_grid, -100.0, 0.5) - - ! Albedo fraction -- read and calculate - call mn_static_read_hires(Atm(1)%npx, Atm(1)%npy, x_refine, Atm(2)%pelist, trim(Moving_nest(child_grid_num)%mn_flag%surface_dir), "facsf", "facsf", mn_static%facsf_grid, parent_tile) - - allocate(mn_static%facwf_grid(lbound(mn_static%facsf_grid,1):ubound(mn_static%facsf_grid,1),lbound(mn_static%facsf_grid,2):ubound(mn_static%facsf_grid,2))) - - ! For land points, set facwf = 1.0 - facsf - ! To match initialization behavior, set any -999s to 0 - do i=lbound(mn_static%facsf_grid,1),ubound(mn_static%facsf_grid,1) - do j=lbound(mn_static%facsf_grid,2),ubound(mn_static%facsf_grid,2) - if (mn_static%facsf_grid(i,j) .lt. -100) then - mn_static%facsf_grid(i,j) = 0 - mn_static%facwf_grid(i,j) = 0 - else - mn_static%facwf_grid(i,j) = 1.0 - mn_static%facsf_grid(i,j) - endif - enddo - enddo + !call validate_geo_coords("LAT", mn_static%fp_ls%geolat_grid, mn_static%nest_ls%geolat_grid, x_refine, ioffset, joffset) + !call validate_geo_coords("LON", mn_static%fp_ls%geolon_grid, mn_static%nest_ls%geolon_grid, x_refine, ioffset, joffset) - ! Additional albedo variables - ! black sky = strong cosz -- direct sunlight - ! white sky = weak cosz -- diffuse light + !! Apply lakes to land mask based on land_frac and soil_type + call mn_apply_lakes(mn_static%parent_ls) + call mn_apply_lakes(mn_static%fp_ls) + call mn_apply_lakes(mn_static%nest_ls) - ! alvsf = visible strong cosz = visible_black_sky_albedo - ! alvwf = visible weak cosz = visible_white_sky_albedo - ! alnsf = near IR strong cosz = near_IR_black_sky_albedo - ! alnwf = near IR weak cosz = near_IR_white_sky_albedo + call mn_static_overwrite_ls_from_nest(mn_static%fp_ls, mn_static%nest_ls, x_refine, ioffset, joffset) - call mn_static_read_hires(Atm(1)%npx, Atm(1)%npy, x_refine, Atm(2)%pelist, trim(Moving_nest(child_grid_num)%mn_flag%surface_dir), "snowfree_albedo", "visible_black_sky_albedo", mn_static%alvsf_grid, parent_tile, time=month) - call mn_static_read_hires(Atm(1)%npx, Atm(1)%npy, x_refine, Atm(2)%pelist, trim(Moving_nest(child_grid_num)%mn_flag%surface_dir), "snowfree_albedo", "visible_white_sky_albedo", mn_static%alvwf_grid, parent_tile, time=month) + ! Initialize the land sea mask (slmsk) in the mn_phys structure + ! Important this is done after adjusting for lakes! + call mn_phys_set_slmsk(Atm, n, mn_static, ioffset, joffset, x_refine) - call mn_static_read_hires(Atm(1)%npx, Atm(1)%npy, x_refine, Atm(2)%pelist, trim(Moving_nest(child_grid_num)%mn_flag%surface_dir), "snowfree_albedo", "near_IR_black_sky_albedo", mn_static%alnsf_grid, parent_tile, time=month) - call mn_static_read_hires(Atm(1)%npx, Atm(1)%npy, x_refine, Atm(2)%pelist, trim(Moving_nest(child_grid_num)%mn_flag%surface_dir), "snowfree_albedo", "near_IR_white_sky_albedo", mn_static%alnwf_grid, parent_tile, time=month) + ! Read in full panel fix data + call mn_static_read_fix(mn_static%fp_fix, Atm(1)%npx, Atm(1)%npy, x_refine, Atm(2)%pelist, trim(Moving_nest(child_grid_num)%mn_flag%surface_dir), parent_tile, month) + ! Read in nest fix data + call mn_static_read_fix(mn_static%nest_fix, Atm(2)%npx, Atm(2)%npy, 1, Atm(2)%pelist, trim(Moving_nest(child_grid_num)%mn_flag%surface_dir) // "/..", static_nest_num, month) - ! Set the -999s to small value of 0.06, matching initialization code in chgres + ! Overwrite fix data from nest initialization + call mn_static_overwrite_fix_from_nest(mn_static%fp_fix, mn_static%nest_fix, x_refine, ioffset, joffset) - call mn_replace_low_values(mn_static%alvsf_grid, -100.0, 0.06) - call mn_replace_low_values(mn_static%alvwf_grid, -100.0, 0.06) - call mn_replace_low_values(mn_static%alnsf_grid, -100.0, 0.06) - call mn_replace_low_values(mn_static%alnwf_grid, -100.0, 0.06) + ! The nest static grids are only used for this step; can safely deallocate them now. + call deallocate_land_mask_grids(mn_static%nest_ls) + call deallocate_fix_grids(mn_static%nest_fix) endif @@ -1007,21 +1218,21 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, select case(Moving_nest(n)%mn_flag%terrain_smoother) case (0) ! High-resolution terrain for entire nest - Atm(n)%phis(isd:ied, jsd:jed) = mn_static%orog_grid((ioffset-1)*x_refine+isd:(ioffset-1)*x_refine+ied, (joffset-1)*y_refine+jsd:(joffset-1)*y_refine+jed) * grav + Atm(n)%phis(isd:ied, jsd:jed) = mn_static%fp_ls%orog_grid((ioffset-1)*x_refine+isd:(ioffset-1)*x_refine+ied, (joffset-1)*y_refine+jsd:(joffset-1)*y_refine+jed) * grav case (1) ! Static nest smoothing algorithm - interpolation of coarse terrain in halo zone and 5 point blending zone of coarse and fine data - call set_blended_terrain(Atm(n), mn_static%parent_orog_grid, mn_static%orog_grid, x_refine, Atm(n)%bd%ng, 5, a_step) + call set_blended_terrain(Atm(n), mn_static%parent_ls%orog_grid, mn_static%fp_ls%orog_grid, x_refine, Atm(n)%bd%ng, 5, a_step) case (2) ! Static nest smoothing algorithm - interpolation of coarse terrain in halo zone and 5 point blending zone of coarse and fine data - call set_blended_terrain(Atm(n), mn_static%parent_orog_grid, mn_static%orog_grid, x_refine, Atm(n)%bd%ng, 10, a_step) + call set_blended_terrain(Atm(n), mn_static%parent_ls%orog_grid, mn_static%fp_ls%orog_grid, x_refine, Atm(n)%bd%ng, 10, a_step) case (4) ! Use coarse terrain; no-op here. ; case (5) ! 5 pt smoother. blend zone of 5 to match static nest - call set_smooth_nest_terrain(Atm(n), mn_static%orog_grid, x_refine, 5, Atm(n)%bd%ng, 5) + call set_smooth_nest_terrain(Atm(n), mn_static%fp_ls%orog_grid, x_refine, 5, Atm(n)%bd%ng, 5) case (9) ! 9 pt smoother. blend zone of 5 to match static nest - call set_smooth_nest_terrain(Atm(n), mn_static%orog_grid, x_refine, 9, Atm(n)%bd%ng, 5) + call set_smooth_nest_terrain(Atm(n), mn_static%fp_ls%orog_grid, x_refine, 9, Atm(n)%bd%ng, 5) case default write (errstring, "(I0)") Moving_nest(n)%mn_flag%terrain_smoother call mpp_error(FATAL,'Invalid terrain_smoother in fv_moving_nest_main '//errstring) @@ -1039,8 +1250,8 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, !real, _ALLOCATABLE :: oro(:,:) _NULL !< land fraction (1: all land; 0: all water) !real, _ALLOCATABLE :: sgh(:,:) _NULL !< Terrain standard deviation - Atm(n)%oro(isc:iec, jsc:jec) = mn_static%land_frac_grid((ioffset-1)*x_refine+isc:(ioffset-1)*x_refine+iec, (joffset-1)*y_refine+jsc:(joffset-1)*y_refine+jec) - Atm(n)%sgh(isc:iec, jsc:jec) = mn_static%orog_std_grid((ioffset-1)*x_refine+isc:(ioffset-1)*x_refine+iec, (joffset-1)*y_refine+jsc:(joffset-1)*y_refine+jec) + Atm(n)%oro(isc:iec, jsc:jec) = mn_static%fp_ls%land_frac_grid((ioffset-1)*x_refine+isc:(ioffset-1)*x_refine+iec, (joffset-1)*y_refine+jsc:(joffset-1)*y_refine+jec) + Atm(n)%sgh(isc:iec, jsc:jec) = mn_static%fp_ls%orog_std_grid((ioffset-1)*x_refine+isc:(ioffset-1)*x_refine+iec, (joffset-1)*y_refine+jsc:(joffset-1)*y_refine+jec) endif call mn_phys_reset_sfc_props(Atm, n, mn_static, Atm_block, IPD_data, ioffset, joffset, x_refine) @@ -1074,15 +1285,39 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, !!===================================================================================== if (use_timers) call mpp_clock_begin (id_movnest7_3) - call mn_prog_apply_temp_variables(Atm, n, child_grid_num, is_fine_pe, npz) - call mn_phys_apply_temp_variables(Atm, Atm_block, IPD_control, IPD_data, n, child_grid_num, is_fine_pe, npz) + if (is_fine_pe) then + ! This just needs to check in the compute domain. While the mn_phys fields extend into the halo, the IPD structure from CCPP only covers the compute domain. + do i=isc,iec + do j=jsc,jec - if (use_timers) call mpp_clock_end (id_movnest7_3) - if (use_timers) call mpp_clock_begin (id_movnest8) + ! Regular physics variables + do k=lbound(Moving_nest(n)%mn_phys%stc,3),ubound(Moving_nest(n)%mn_phys%stc,3) + if (Moving_nest(n)%mn_phys%stc(i,j,k) .lt. 0.0 .or. Moving_nest(n)%mn_phys%stc(i,j,k) .gt. maxSkinTempK ) then + print '("[INFO] WDR PHYSICS reset soil temp values npe=",I0," i=",I0," j=",I0," stc=",E12.5)', mpp_pe(), i, j, Moving_nest(n)%mn_phys%stc(i,j,k) + Moving_nest(n)%mn_phys%stc(i,j,k) = 298.0 + endif + if (Moving_nest(n)%mn_phys%smc(i,j,k) .lt. 0.0 .or. Moving_nest(n)%mn_phys%smc(i,j,k) .gt. 1000.0 ) then + print '("[INFO] WDR PHYSICS reset soil moisture values npe=",I0," i=",I0," j=",I0," smc=",E12.5)', mpp_pe(), i, j, Moving_nest(n)%mn_phys%smc(i,j,k) + Moving_nest(n)%mn_phys%smc(i,j,k) = 0.3 + endif + if (Moving_nest(n)%mn_phys%slc(i,j,k) .lt. 0.0 .or. Moving_nest(n)%mn_phys%slc(i,j,k) .gt. 1000.0 ) then + print '("[INFO] WDR PHYSICS reset soil liquid values npe=",I0," i=",I0," j=",I0," slc=",E12.5)', mpp_pe(), i, j, Moving_nest(n)%mn_phys%slc(i,j,k) + Moving_nest(n)%mn_phys%slc(i,j,k) = 0.3 + endif + enddo - !!============================================================================ - !! Step 8 -- Dump to netCDF - !!============================================================================ + + if (Moving_nest(n)%mn_phys%canopy(i,j) .lt. 0.0 .or. Moving_nest(n)%mn_phys%canopy(i,j) .gt. 100.0 ) then + print '("[INFO] WDR PHYSICS reset canopy water values npe=",I0," i=",I0," j=",I0," canopy=",E12.5)', mpp_pe(), i, j, Moving_nest(n)%mn_phys%canopy(i,j) + Moving_nest(n)%mn_phys%canopy(i,j) = 0.0 ! Zero out if no other information + endif + if (Moving_nest(n)%mn_phys%vegfrac(i,j) .lt. 0.0 .or. Moving_nest(n)%mn_phys%vegfrac(i,j) .gt. 100.0 ) then + print '("[INFO] WDR PHYSICS reset vegfrac values npe=",I0," i=",I0," j=",I0," vegfrac=",E12.5)', mpp_pe(), i, j, Moving_nest(n)%mn_phys%vegfrac(i,j) + Moving_nest(n)%mn_phys%vegfrac(i,j) = 0.5 ! Default to half. Confirmed that values are fractions. + endif + enddo ! do j + enddo ! do i + endif ! if is_fine_pe if (is_fine_pe) then @@ -1106,10 +1341,21 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, if (Moving_nest(n)%mn_phys%albdifvis_lnd(i,j) .lt. 0.0) Moving_nest(n)%mn_phys%albdifvis_lnd(i,j) = 0.5 if (Moving_nest(n)%mn_phys%albdifnir_lnd(i,j) .lt. 0.0) Moving_nest(n)%mn_phys%albdifnir_lnd(i,j) = 0.5 + enddo enddo endif + call mn_prog_apply_temp_variables(Atm, n, child_grid_num, is_fine_pe, npz) + call mn_phys_apply_temp_variables(Atm, Atm_block, IPD_control, IPD_data, n, child_grid_num, is_fine_pe, npz, a_step) + + if (use_timers) call mpp_clock_end (id_movnest7_3) + if (use_timers) call mpp_clock_begin (id_movnest8) + + !!============================================================================ + !! Step 8 -- Dump to netCDF + !!============================================================================ + output_step = output_step + 1 if (debug_sync) call mpp_sync(full_pelist) ! Used to make debugging easier. Can be removed. @@ -1148,20 +1394,4 @@ subroutine fv_moving_nest_exec(Atm, Atm_block, IPD_control, IPD_data, delta_i_c, end subroutine fv_moving_nest_exec - !>@brief The subroutine 'mn_replace_low_values' replaces low values with a default value. - subroutine mn_replace_low_values(data_grid, low_value, new_value) - real, _ALLOCATABLE, intent(inout) :: data_grid(:,:) !< 2D grid of data - real, intent(in) :: low_value !< Low value to check for; e.g. negative or fill value - real, intent(in) :: new_value !< Value to replace low value with - - integer :: i, j - - do i=lbound(data_grid,1),ubound(data_grid,1) - do j=lbound(data_grid,2),ubound(data_grid,2) - if (data_grid(i,j) .le. low_value) data_grid(i,j) = new_value - enddo - enddo - end subroutine mn_replace_low_values - end module fv_moving_nest_main_mod - diff --git a/moving_nest/fv_moving_nest_physics.F90 b/moving_nest/fv_moving_nest_physics.F90 index 1219617a6..f49b024f6 100644 --- a/moving_nest/fv_moving_nest_physics.F90 +++ b/moving_nest/fv_moving_nest_physics.F90 @@ -123,6 +123,29 @@ module fv_moving_nest_physics_mod contains + + subroutine mn_phys_set_slmsk(Atm, n, mn_static, ioffset, joffset, refine) + type(fv_atmos_type), intent(inout),allocatable :: Atm(:) !< Array of atmospheric data + integer, intent(in) :: n !< Current grid number + type(mn_surface_grids), intent(in) :: mn_static !< Static surface data + integer, intent(in) :: ioffset, joffset !< Current nest offset in i,j direction + integer, intent(in) :: refine !< Nest refinement ratio + + integer :: i_pe, j_pe, i_idx, j_idx + + !print '("[INFO] MASK inside mn_phys_set_slmsk npe=",I0)', mpp_pe() + ! Setup local land sea mask grid for masked interpolations + do i_pe = Atm(n)%bd%isd, Atm(n)%bd%ied + do j_pe = Atm(n)%bd%jsd, Atm(n)%bd%jed + i_idx = (ioffset-1)*refine + i_pe + j_idx = (joffset-1)*refine + j_pe + + Moving_nest(n)%mn_phys%slmsk(i_pe, j_pe) = mn_static%fp_ls%ls_mask_grid(i_idx, j_idx) + enddo + enddo + end subroutine mn_phys_set_slmsk + + !>@brief The subroutine 'mn_phys_reset_sfc_props' sets the static surface parameters from the high-resolution input file data !>@details This subroutine relies on earlier code reading the data from files into the mn_static data structure !! This subroutine does not yet handle ice points or frac_grid - fractional landfrac/oceanfrac values @@ -139,15 +162,8 @@ subroutine mn_phys_reset_sfc_props(Atm, n, mn_static, Atm_block, IPD_data, ioffs integer :: nb, blen, ix, i_pe, j_pe, i_idx, j_idx real(kind=kind_phys) :: phys_oro - ! Setup local land sea mask grid for masked interpolations - do i_pe = Atm(n)%bd%isd, Atm(n)%bd%ied - do j_pe = Atm(n)%bd%jsd, Atm(n)%bd%jed - i_idx = (ioffset-1)*refine + i_pe - j_idx = (joffset-1)*refine + j_pe - - Moving_nest(n)%mn_phys%slmsk(i_pe, j_pe) = mn_static%ls_mask_grid(i_idx, j_idx) - enddo - enddo + !print '("[INFO] MASK inside mn_phys_reset_sfc_props npe=",I0)', mpp_pe() + call mn_phys_set_slmsk(Atm, n, mn_static, ioffset, joffset, refine) ! Reset the variables from the fix_sfc files do nb = 1,Atm_block%nblks @@ -160,12 +176,12 @@ subroutine mn_phys_reset_sfc_props(Atm, n, mn_static, Atm_block, IPD_data, ioffs j_idx = (joffset-1)*refine + j_pe ! Reset the land sea mask from the hires parent data - IPD_data(nb)%Sfcprop%slmsk(ix) = mn_static%ls_mask_grid(i_idx, j_idx) + IPD_data(nb)%Sfcprop%slmsk(ix) = mn_static%fp_ls%ls_mask_grid(i_idx, j_idx) ! IFD values are 0 for land, and 1 for oceans/lakes -- reverse of the land sea mask ! Land Sea Mask has values of 0 for oceans/lakes, 1 for land, 2 for sea ice ! TODO figure out what ifd should be for sea ice - if (mn_static%ls_mask_grid(i_idx, j_idx) .eq. 1 ) then + if (mn_static%fp_ls%ls_mask_grid(i_idx, j_idx) .eq. 1 ) then if (move_nsst) IPD_data(nb)%Sfcprop%ifd(ix) = 0 ! Land IPD_data(nb)%Sfcprop%oceanfrac(ix) = 0 ! Land -- TODO permit fractions IPD_data(nb)%Sfcprop%landfrac(ix) = 1 ! Land -- TODO permit fractions @@ -175,16 +191,16 @@ subroutine mn_phys_reset_sfc_props(Atm, n, mn_static, Atm_block, IPD_data, ioffs IPD_data(nb)%Sfcprop%landfrac(ix) = 0 ! Ocean -- TODO permit fractions endif - IPD_data(nb)%Sfcprop%tg3(ix) = mn_static%deep_soil_temp_grid(i_idx, j_idx) + IPD_data(nb)%Sfcprop%tg3(ix) = mn_static%fp_fix%deep_soil_temp_grid(i_idx, j_idx) - ! Follow logic from FV3/io/fv3atm_sfc_io.F90 + ! Follow logic from FV3/io/FV3GFS_io.F90 line 1187 ! TODO this will need to be more complicated if we support frac_grid !if (nint(mn_static%soil_type_grid(i_idx, j_idx)) == 14 .or. int(mn_static%soil_type_grid(i_idx, j_idx)+0.5) <= 0) then !if (nint(mn_static%soil_type_grid(i_idx, j_idx)) == 14 .or. !if ( (mn_static%ls_mask_grid(i_idx, j_idx) .eq. 1 .and. nint(mn_static%land_frac_grid(i_idx, j_idx)) == 0) .or. & ! mn_static%soil_type_grid(i_idx, j_idx) < 0.5) then - if (mn_static%ls_mask_grid(i_idx, j_idx) .eq. 1 .and. nint(mn_static%land_frac_grid(i_idx, j_idx)) == 0 ) then + if (mn_static%fp_ls%ls_mask_grid(i_idx, j_idx) .eq. 1 .and. nint(mn_static%fp_ls%land_frac_grid(i_idx, j_idx)) == 0 ) then ! Water soil type == lake, etc. -- override the other variables and make this water !!print '("mn_phys_reset_sfc_props LAKE SOIL npe=",I0," x,y=",I0,",",I0," lat=",F10.3," lon=",F10.3)', mpp_pe(), i_idx, j_idx, IPD_data(nb)%Grid%xlat_d(ix), IPD_data(nb)%Grid%xlon_d(ix)-360.0 @@ -195,21 +211,21 @@ subroutine mn_phys_reset_sfc_props(Atm, n, mn_static, Atm_block, IPD_data, ioffs IPD_data(nb)%Sfcprop%stype(ix) = 0 IPD_data(nb)%Sfcprop%slmsk(ix) = 0 else - IPD_data(nb)%Sfcprop%stype(ix) = nint(mn_static%soil_type_grid(i_idx, j_idx)) + IPD_data(nb)%Sfcprop%stype(ix) = nint(mn_static%fp_ls%soil_type_grid(i_idx, j_idx)) endif !IPD_data(nb)%Sfcprop%vfrac(ix) = mn_static%veg_frac_grid(i_idx, j_idx) - IPD_data(nb)%Sfcprop%vtype(ix) = nint(mn_static%veg_type_grid(i_idx, j_idx)) - IPD_data(nb)%Sfcprop%slope(ix) = nint(mn_static%slope_type_grid(i_idx, j_idx)) - IPD_data(nb)%Sfcprop%snoalb(ix) = mn_static%max_snow_alb_grid(i_idx, j_idx) + IPD_data(nb)%Sfcprop%vtype(ix) = nint(mn_static%fp_fix%veg_type_grid(i_idx, j_idx)) + IPD_data(nb)%Sfcprop%slope(ix) = nint(mn_static%fp_fix%slope_type_grid(i_idx, j_idx)) + IPD_data(nb)%Sfcprop%snoalb(ix) = mn_static%fp_fix%max_snow_alb_grid(i_idx, j_idx) - IPD_data(nb)%Sfcprop%facsf(ix) = mn_static%facsf_grid(i_idx, j_idx) - IPD_data(nb)%Sfcprop%facwf(ix) = mn_static%facwf_grid(i_idx, j_idx) + IPD_data(nb)%Sfcprop%facsf(ix) = mn_static%fp_fix%facsf_grid(i_idx, j_idx) + IPD_data(nb)%Sfcprop%facwf(ix) = mn_static%fp_fix%facwf_grid(i_idx, j_idx) - IPD_data(nb)%Sfcprop%alvsf(ix) = mn_static%alvsf_grid(i_idx, j_idx) - IPD_data(nb)%Sfcprop%alvwf(ix) = mn_static%alvwf_grid(i_idx, j_idx) - IPD_data(nb)%Sfcprop%alnsf(ix) = mn_static%alnsf_grid(i_idx, j_idx) - IPD_data(nb)%Sfcprop%alnwf(ix) = mn_static%alnwf_grid(i_idx, j_idx) + IPD_data(nb)%Sfcprop%alvsf(ix) = mn_static%fp_fix%alvsf_grid(i_idx, j_idx) + IPD_data(nb)%Sfcprop%alvwf(ix) = mn_static%fp_fix%alvwf_grid(i_idx, j_idx) + IPD_data(nb)%Sfcprop%alnsf(ix) = mn_static%fp_fix%alnsf_grid(i_idx, j_idx) + IPD_data(nb)%Sfcprop%alnwf(ix) = mn_static%fp_fix%alnwf_grid(i_idx, j_idx) ! Reset the orography in the physics arrays, using the smoothed values from above phys_oro = Atm(n)%phis(i_pe, j_pe) / grav @@ -313,6 +329,7 @@ subroutine mn_phys_fill_temp_variables(Atm, Atm_block, IPD_Control, IPD_Data, n, integer :: nb, blen, i, j, k, ix, nv type(fv_moving_nest_physics_type), pointer :: mn_phys + integer :: err_field = 0 this_pe = mpp_pe() @@ -438,7 +455,7 @@ end subroutine mn_phys_fill_temp_variables !>@brief The subroutine 'mn_phys_apply_temp_variables' copies moved 2D data back into 1D physics arryas for nest motion !>@details This subroutine fills the 1D physics arrays from the mn_phys structure on the Atm object !! Note that ice variables are not yet handled. - subroutine mn_phys_apply_temp_variables(Atm, Atm_block, IPD_Control, IPD_Data, n, child_grid_num, is_fine_pe, npz) + subroutine mn_phys_apply_temp_variables(Atm, Atm_block, IPD_Control, IPD_Data, n, child_grid_num, is_fine_pe, npz, a_step) type(fv_atmos_type), allocatable, target, intent(inout) :: Atm(:) !< Array of atmospheric data type (block_control_type), intent(in) :: Atm_block !< Physics block layout type(IPD_control_type), intent(in) :: IPD_Control !< Physics metadata @@ -446,6 +463,7 @@ subroutine mn_phys_apply_temp_variables(Atm, Atm_block, IPD_Control, IPD_Data, n integer, intent(in) :: n, child_grid_num !< Current grid number, child grid number logical, intent(in) :: is_fine_pe !< Is this a nest PE? integer, intent(in) :: npz !< Number of vertical levels + integer, intent(in) :: a_step !< Timestep for logging integer :: is, ie, js, je integer :: this_pe @@ -474,7 +492,13 @@ subroutine mn_phys_apply_temp_variables(Atm, Atm_block, IPD_Control, IPD_Data, n if (move_physics) then ! Surface properties + !print '("[INFO] WDR smc set npe=",I0," smc(",I0,",",I0,",",I0,")=",E18.10," slmsk=",I0)', this_pe, i, j, 1, mn_phys%smc(i,j,1), int(mn_phys%slmsk(i,j)) do k = 1, IPD_Control%lsoil + !if ( int(mn_phys%slmsk(i,j)) .eq. 1 .and. mn_phys%smc(i,j,k) .ge. 100) then + ! IPD_Data(nb)%Sfcprop%smc(ix,k) = 0.3 + !else + ! IPD_Data(nb)%Sfcprop%smc(ix,k) = mn_phys%smc(i,j,k) + !endif IPD_Data(nb)%Sfcprop%smc(ix,k) = mn_phys%smc(i,j,k) IPD_Data(nb)%Sfcprop%stc(ix,k) = mn_phys%stc(i,j,k) IPD_Data(nb)%Sfcprop%slc(ix,k) = mn_phys%slc(i,j,k) @@ -610,26 +634,6 @@ subroutine mn_phys_apply_temp_variables(Atm, Atm_block, IPD_Control, IPD_Data, n IPD_Data(nb)%Sfcprop%dt_cool(ix) = mn_phys%dt_cool(i,j) IPD_Data(nb)%Sfcprop%qrain(ix) = mn_phys%qrain(i,j) endif - - ! Check if stype and vtype are properly set for land points. Set to reasonable values if they have fill values. - if ( (int(IPD_data(nb)%Sfcprop%slmsk(ix)) .eq. 1) ) then - - if (IPD_data(nb)%Sfcprop%vtype(ix) .lt. 0.5) then - IPD_data(nb)%Sfcprop%vtype(ix) = 7 ! Force to grassland - endif - - if (IPD_data(nb)%Sfcprop%stype(ix) .lt. 0.5) then - IPD_data(nb)%Sfcprop%stype(ix) = 3 ! Force to sandy loam - endif - - if (IPD_data(nb)%Sfcprop%vtype_save(ix) .lt. 0.5) then - IPD_data(nb)%Sfcprop%vtype_save(ix) = 7 ! Force to grassland - endif - if (IPD_data(nb)%Sfcprop%stype_save(ix) .lt. 0.5) then - IPD_data(nb)%Sfcprop%stype_save(ix) = 3 ! Force to sandy loam - endif - - endif enddo enddo endif @@ -654,6 +658,8 @@ subroutine mn_phys_fill_nest_halos_from_parent(Atm, IPD_Control, IPD_Data, mn_st integer :: x_refine, y_refine type(fv_moving_nest_physics_type), pointer :: mn_phys + integer, parameter :: M_WATER = 0, M_LAND = 1, M_SEAICE = 2 + interp_type = 1 ! cell-centered A-grid interp_type_u = 4 ! D-grid interp_type_v = 4 ! D-grid @@ -676,18 +682,22 @@ subroutine mn_phys_fill_nest_halos_from_parent(Atm, IPD_Control, IPD_Data, mn_st is_fine_pe, nest_domain, position) if (move_physics) then - call fill_nest_halos_from_parent("smc", mn_phys%smc, interp_type, Atm(child_grid_num)%neststruct%wt_h, & + ! Default - Arbitrary value 0.3 + call fill_nest_halos_from_parent_masked("smc", mn_phys%smc, interp_type_lmask, Atm(child_grid_num)%neststruct%wt_h, & Atm(child_grid_num)%neststruct%ind_h, & x_refine, y_refine, & - is_fine_pe, nest_domain, position, IPD_Control%lsoil) - call fill_nest_halos_from_parent("stc", mn_phys%stc, interp_type, Atm(child_grid_num)%neststruct%wt_h, & + is_fine_pe, nest_domain, position, 1, IPD_Control%lsoil, mn_phys%slmsk, mn_static%parent_ls%ls_mask_grid, M_LAND, 0.3D0) + ! Defaults - use surface temperature to set soil temperature at each level + call fill_nest_halos_from_parent_masked("stc", mn_phys%stc, interp_type_lmask, Atm(child_grid_num)%neststruct%wt_h, & Atm(child_grid_num)%neststruct%ind_h, & x_refine, y_refine, & - is_fine_pe, nest_domain, position, IPD_Control%lsoil) - call fill_nest_halos_from_parent("slc", mn_phys%slc, interp_type, Atm(child_grid_num)%neststruct%wt_h, & + is_fine_pe, nest_domain, position, 1, IPD_Control%lsoil, mn_phys%slmsk, mn_static%parent_ls%ls_mask_grid, M_LAND, mn_phys%ts) + ! Default - Arbitrary value 0.3 + call fill_nest_halos_from_parent_masked("slc", mn_phys%slc, interp_type_lmask, Atm(child_grid_num)%neststruct%wt_h, & Atm(child_grid_num)%neststruct%ind_h, & x_refine, y_refine, & - is_fine_pe, nest_domain, position, IPD_Control%lsoil) + is_fine_pe, nest_domain, position, 1, IPD_Control%lsoil, mn_phys%slmsk, mn_static%parent_ls%ls_mask_grid, M_LAND, 0.3D0) + call fill_nest_halos_from_parent("phy_f2d", mn_phys%phy_f2d, interp_type, Atm(child_grid_num)%neststruct%wt_h, & Atm(child_grid_num)%neststruct%ind_h, & @@ -707,21 +717,22 @@ subroutine mn_phys_fill_nest_halos_from_parent(Atm, IPD_Control, IPD_Data, mn_st ! is_fine_pe, nest_domain, position) ! sea/land mask array (sea:0,land:1,sea-ice:2) + !integer, parameter :: M_WATER = 0, M_LAND = 1, M_SEAICE = 2 call fill_nest_halos_from_parent_masked("emis_lnd", mn_phys%emis_lnd, interp_type_lmask, Atm(child_grid_num)%neststruct%wt_h, & Atm(child_grid_num)%neststruct%ind_h, & x_refine, y_refine, & - is_fine_pe, nest_domain, position, mn_phys%slmsk, 1, 0.5D0) + is_fine_pe, nest_domain, position, mn_phys%slmsk, mn_static%parent_ls%ls_mask_grid, M_LAND, 0.5D0) call fill_nest_halos_from_parent_masked("emis_ice", mn_phys%emis_ice, interp_type_lmask, Atm(child_grid_num)%neststruct%wt_h, & Atm(child_grid_num)%neststruct%ind_h, & x_refine, y_refine, & - is_fine_pe, nest_domain, position, mn_phys%slmsk, 2, 0.5D0) + is_fine_pe, nest_domain, position, mn_phys%slmsk, mn_static%parent_ls%ls_mask_grid, M_SEAICE, 0.5D0) call fill_nest_halos_from_parent_masked("emis_wat", mn_phys%emis_wat, interp_type_lmask, Atm(child_grid_num)%neststruct%wt_h, & Atm(child_grid_num)%neststruct%ind_h, & x_refine, y_refine, & - is_fine_pe, nest_domain, position, mn_phys%slmsk, 0, 0.5D0) + is_fine_pe, nest_domain, position, mn_phys%slmsk, mn_static%parent_ls%ls_mask_grid, M_WATER, 0.5D0) !call fill_nest_halos_from_parent("sfalb_lnd_bck", mn_phys%sfalb_lnd_bck, interp_type, Atm(child_grid_num)%neststruct%wt_h, & ! Atm(child_grid_num)%neststruct%ind_h, & @@ -770,14 +781,23 @@ subroutine mn_phys_fill_nest_halos_from_parent(Atm, IPD_Control, IPD_Data, mn_st x_refine, y_refine, & is_fine_pe, nest_domain, position) - call fill_nest_halos_from_parent("canopy", mn_phys%canopy, interp_type, Atm(child_grid_num)%neststruct%wt_h, & - Atm(child_grid_num)%neststruct%ind_h, & - x_refine, y_refine, & - is_fine_pe, nest_domain, position) - call fill_nest_halos_from_parent("vegfrac", mn_phys%vegfrac, interp_type, Atm(child_grid_num)%neststruct%wt_h, & - Atm(child_grid_num)%neststruct%ind_h, & - x_refine, y_refine, & - is_fine_pe, nest_domain, position) +! call fill_nest_halos_from_parent("canopy", mn_phys%canopy, interp_type, Atm(child_grid_num)%neststruct%wt_h, & +! Atm(child_grid_num)%neststruct%ind_h, & +! x_refine, y_refine, & +! is_fine_pe, nest_domain, position) +! call fill_nest_halos_from_parent("vegfrac", mn_phys%vegfrac, interp_type, Atm(child_grid_num)%neststruct%wt_h, & +! Atm(child_grid_num)%neststruct%ind_h, & +! x_refine, y_refine, & +! is_fine_pe, nest_domain, position) + call fill_nest_halos_from_parent_masked("canopy", mn_phys%canopy, interp_type_lmask, Atm(child_grid_num)%neststruct%wt_h, & + Atm(child_grid_num)%neststruct%ind_h, x_refine, y_refine, & + is_fine_pe, nest_domain, position, mn_phys%slmsk, mn_static%parent_ls%ls_mask_grid, M_LAND, 0.0D0) + call fill_nest_halos_from_parent_masked("vegfrac", mn_phys%vegfrac, interp_type_lmask, Atm(child_grid_num)%neststruct%wt_h, & + Atm(child_grid_num)%neststruct%ind_h, x_refine, y_refine, & + is_fine_pe, nest_domain, position, mn_phys%slmsk, mn_static%parent_ls%ls_mask_grid, M_LAND, 0.50D0) + + + call fill_nest_halos_from_parent("uustar", mn_phys%uustar, interp_type, Atm(child_grid_num)%neststruct%wt_h, & Atm(child_grid_num)%neststruct%ind_h, & x_refine, y_refine, & @@ -798,15 +818,15 @@ subroutine mn_phys_fill_nest_halos_from_parent(Atm, IPD_Control, IPD_Data, mn_st call fill_nest_halos_from_parent_masked("zorll", mn_phys%zorll, interp_type_lmask, Atm(child_grid_num)%neststruct%wt_h, & Atm(child_grid_num)%neststruct%ind_h, & x_refine, y_refine, & - is_fine_pe, nest_domain, position, mn_phys%slmsk, 1, 86.0D0) + is_fine_pe, nest_domain, position, mn_phys%slmsk, mn_static%parent_ls%ls_mask_grid, 1, 86.0D0) call fill_nest_halos_from_parent_masked("zorlwav", mn_phys%zorlwav, interp_type_lmask, Atm(child_grid_num)%neststruct%wt_h, & Atm(child_grid_num)%neststruct%ind_h, & x_refine, y_refine, & - is_fine_pe, nest_domain, position, mn_phys%slmsk, 0, 77.0D0) + is_fine_pe, nest_domain, position, mn_phys%slmsk, mn_static%parent_ls%ls_mask_grid, 0, 77.0D0) call fill_nest_halos_from_parent_masked("zorlw", mn_phys%zorlw, interp_type_lmask, Atm(child_grid_num)%neststruct%wt_h, & Atm(child_grid_num)%neststruct%ind_h, & x_refine, y_refine, & - is_fine_pe, nest_domain, position, mn_phys%slmsk, 0, 78.0D0) + is_fine_pe, nest_domain, position, mn_phys%slmsk, mn_static%parent_ls%ls_mask_grid, 0, 78.0D0) call fill_nest_halos_from_parent("tsfco", mn_phys%tsfco, interp_type, Atm(child_grid_num)%neststruct%wt_h, & Atm(child_grid_num)%neststruct%ind_h, & @@ -824,19 +844,19 @@ subroutine mn_phys_fill_nest_halos_from_parent(Atm, IPD_Control, IPD_Data, mn_st call fill_nest_halos_from_parent_masked("albdirvis_lnd", mn_phys%albdirvis_lnd, interp_type_lmask, Atm(child_grid_num)%neststruct%wt_h, & Atm(child_grid_num)%neststruct%ind_h, & x_refine, y_refine, & - is_fine_pe, nest_domain, position, mn_phys%slmsk, 1, 0.5D0) + is_fine_pe, nest_domain, position, mn_phys%slmsk, mn_static%parent_ls%ls_mask_grid, M_LAND, 0.5D0) call fill_nest_halos_from_parent_masked("albdirnir_lnd", mn_phys%albdirnir_lnd, interp_type_lmask, Atm(child_grid_num)%neststruct%wt_h, & Atm(child_grid_num)%neststruct%ind_h, & x_refine, y_refine, & - is_fine_pe, nest_domain, position, mn_phys%slmsk, 1, 0.5D0) + is_fine_pe, nest_domain, position, mn_phys%slmsk, mn_static%parent_ls%ls_mask_grid, M_LAND, 0.5D0) call fill_nest_halos_from_parent_masked("albdifvis_lnd", mn_phys%albdifvis_lnd, interp_type_lmask, Atm(child_grid_num)%neststruct%wt_h, & Atm(child_grid_num)%neststruct%ind_h, & x_refine, y_refine, & - is_fine_pe, nest_domain, position, mn_phys%slmsk, 1, 0.5D0) + is_fine_pe, nest_domain, position, mn_phys%slmsk, mn_static%parent_ls%ls_mask_grid, M_LAND, 0.5D0) call fill_nest_halos_from_parent_masked("albdifnir_lnd", mn_phys%albdifnir_lnd, interp_type_lmask, Atm(child_grid_num)%neststruct%wt_h, & Atm(child_grid_num)%neststruct%ind_h, & x_refine, y_refine, & - is_fine_pe, nest_domain, position, mn_phys%slmsk, 1, 0.5D0) + is_fine_pe, nest_domain, position, mn_phys%slmsk, mn_static%parent_ls%ls_mask_grid, M_LAND, 0.5D0) @@ -1211,6 +1231,7 @@ subroutine mn_phys_dump_to_netcdf(Atm, Atm_block, IPD_Control, IPD_Data, time_va real, allocatable :: lakefrac_pr_local (:,:) !< lake fraction real, allocatable :: landfrac_pr_local (:,:) !< land fraction real, allocatable :: emis_lnd_pr_local (:,:) !< emissivity land + real, allocatable :: snowxy_pr_local (:,:) !< number of snow layers this_pe = mpp_pe() @@ -1302,6 +1323,9 @@ subroutine mn_phys_dump_to_netcdf(Atm, Atm_block, IPD_Control, IPD_Data, time_va xv_pr_local = +99999.9 ifd_pr_local = +99999.9 endif + if (move_nsst) then + snowxy_pr_local = +99999.9 + endif do nb = 1,Atm_block%nblks blen = Atm_block%blksz(nb) diff --git a/moving_nest/fv_moving_nest_types.F90 b/moving_nest/fv_moving_nest_types.F90 index 45bd12504..5d7a00c04 100644 --- a/moving_nest/fv_moving_nest_types.F90 +++ b/moving_nest/fv_moving_nest_types.F90 @@ -76,26 +76,30 @@ module fv_moving_nest_types_mod real, _ALLOCATABLE :: delz(:,:,:) _NULL !< layer thickness (meters) end type fv_moving_nest_prog_type - ! TODO deallocate these at end of model run. They are only allocated once, at first nest move, inside mn_static_read_hires(). - ! Note these are only 32 bits for now; matching the precision of the input netCDF files - ! though the model generally handles physics variables with 64 bit precision - type mn_surface_grids + + type mn_land_mask_grids real, allocatable :: orog_grid(:,:) _NULL ! orography -- raw or filtered depending on namelist option, in meters real, allocatable :: orog_std_grid(:,:) _NULL ! terrain standard deviation for gravity wave drag, in meters (?) real, allocatable :: ls_mask_grid(:,:) _NULL ! land sea mask -- 0 for ocean/lakes, 1, for land. Perhaps 2 for sea ice. - real, allocatable :: land_frac_grid(:,:) _NULL ! Continuous land fraction - 0.0 ocean, 0.5 half of each, 1.0 all land + real, allocatable :: soil_type_grid(:,:) _NULL ! STATSGO soil type + ! Land frac needs to be kind_phys because CCPP defines it that way. Can have rounding mismatches around 0.5 if types don't match. + real(kind=kind_phys), allocatable :: land_frac_grid(:,:) _NULL ! Continuous land fraction - 0.0 ocean, 0.5 half of each, 1.0 all land - real, allocatable :: parent_orog_grid(:,:) _NULL ! parent orography -- only used for terrain_smoother=1. ! raw or filtered depending on namelist option,in meters + real(kind=kind_phys), allocatable :: geolat_grid(:,:) _NULL + real(kind=kind_phys), allocatable :: geolon_grid(:,:) _NULL + end type mn_land_mask_grids + + type mn_fix_grids ! Soil variables real, allocatable :: deep_soil_temp_grid(:,:) _NULL ! deep soil temperature at 5m, in degrees K - real, allocatable :: soil_type_grid(:,:) _NULL ! STATSGO soil type ! Vegetation variables real, allocatable :: veg_frac_grid(:,:) _NULL ! vegetation fraction real, allocatable :: veg_type_grid(:,:) _NULL ! IGBP vegetation type - real, allocatable :: veg_greenness_grid(:,:) _NULL ! NESDIS vegetation greenness; netCDF file has monthly values + ! TODO do we need veg_greenness? + !real, allocatable :: veg_greenness_grid(:,:) _NULL ! NESDIS vegetation greenness; netCDF file has monthly values ! Orography variables real, allocatable :: slope_type_grid(:,:) _NULL ! legacy 1 degree GFS slope type @@ -118,6 +122,22 @@ module fv_moving_nest_types_mod real, allocatable :: alvwf_grid(:,:) _NULL ! Visible white sky albedo; netCDF file has monthly values real, allocatable :: alnsf_grid(:,:) _NULL ! Near IR black sky albedo; netCDF file has monthly values real, allocatable :: alnwf_grid(:,:) _NULL ! Near IR white sky albedo; netCDF file has monthly values + end type mn_fix_grids + + + + ! TODO deallocate these at end of model run. They are only allocated once, at first nest move, inside mn_static_read_hires(). + ! Note these are only 32 bits for now; matching the precision of the input netCDF files + ! though the model generally handles physics variables with 64 bit precision + type mn_surface_grids + + type(mn_land_mask_grids) :: parent_ls + type(mn_land_mask_grids) :: fp_ls + type(mn_land_mask_grids) :: nest_ls + + ! type(mn_fix_grids) :: parent_fix ! Not needed at present + type(mn_fix_grids) :: fp_fix + type(mn_fix_grids) :: nest_fix end type mn_surface_grids @@ -243,6 +263,12 @@ module fv_moving_nest_types_mod type(fv_moving_nest_type), _ALLOCATABLE, target :: Moving_nest(:) + interface mn_overwrite_with_nest_init_values + module procedure mn_overwrite_with_nest_init_values_r4 + module procedure mn_overwrite_with_nest_init_values_r8 + end interface mn_overwrite_with_nest_init_values + + contains subroutine fv_moving_nest_init(Atm, this_grid) @@ -332,6 +358,169 @@ subroutine deallocate_fv_moving_nest(n) end subroutine deallocate_fv_moving_nest + subroutine mn_apply_lakes(land_mask_grids) + type(mn_land_mask_grids), intent(inout) :: land_mask_grids + + integer :: i_idx, j_idx + + ! Alter hires full panel ls_mask_grid to set lakes to water(sea) values + do i_idx = lbound(land_mask_grids%ls_mask_grid,1), ubound(land_mask_grids%ls_mask_grid,1) + do j_idx = lbound(land_mask_grids%ls_mask_grid,1), ubound(land_mask_grids%ls_mask_grid,2) + !if (land_mask_grids%ls_mask_grid(i_idx, j_idx) .eq. 1 .and. nint(land_mask_grids%land_frac_grid(i_idx, j_idx)) == 0 ) then + !!if (land_mask_grids%ls_mask_grid(i_idx, j_idx) .eq. 1 .and. land_mask_grids%land_frac_grid(i_idx, j_idx) .lt. 0.999 ) then + + ! Use epsilon of 1.0e-6 on land_frac_grid, based on CCPP code in physics/physics/gcycle.F90 + ! Fixes a bug where land mask changes with first nest move if land_frac_grid = 0.5000 + ! TODO test wrapping these reals with int() or nint() + if (land_mask_grids%ls_mask_grid(i_idx, j_idx) .eq. 1 .and. nint(land_mask_grids%land_frac_grid(i_idx, j_idx)-1.0e-6_kind_phys) .eq. 0 ) then + land_mask_grids%ls_mask_grid(i_idx, j_idx) = 0 + endif + ! Soil type adjustments from io/fv3atm_sfc_io.F90 + if (land_mask_grids%ls_mask_grid(i_idx, j_idx) .eq. 1 .and. int(land_mask_grids%soil_type_grid(i_idx, j_idx)) .eq. 14 ) then + land_mask_grids%ls_mask_grid(i_idx, j_idx) = 0 + endif + if (land_mask_grids%ls_mask_grid(i_idx, j_idx) .eq. 1 .and. land_mask_grids%soil_type_grid(i_idx, j_idx) .lt. 0.8 ) then + land_mask_grids%ls_mask_grid(i_idx, j_idx) = 0 + endif + enddo + enddo + + end subroutine mn_apply_lakes + + subroutine mn_overwrite_with_nest_init_values_r8(tag, var_grid, nest_var_grid, refine, ioffset, joffset) + character(len=*) :: tag + real*8, allocatable, intent(inout) :: var_grid(:,:) + real*8, allocatable, intent(in) :: nest_var_grid(:,:) + + integer, intent(in) :: refine, ioffset, joffset + integer :: i,j, this_pe + + ! this_pe = mpp_pe() + + do i = lbound(nest_var_grid,1), ubound(nest_var_grid,1) + do j = lbound(nest_var_grid,2), ubound(nest_var_grid,2) + var_grid((ioffset-1)*refine+i, (joffset-1)*refine+j) = nest_var_grid(i,j) + enddo + enddo + + end subroutine mn_overwrite_with_nest_init_values_r8 + + subroutine mn_overwrite_with_nest_init_values_r4(tag, var_grid, nest_var_grid, refine, ioffset, joffset) + character(len=*) :: tag + real*4, allocatable, intent(inout) :: var_grid(:,:) + real*4, allocatable, intent(in) :: nest_var_grid(:,:) + + integer, intent(in) :: refine, ioffset, joffset + integer :: i,j, this_pe + + ! this_pe = mpp_pe() + + do i = lbound(nest_var_grid,1), ubound(nest_var_grid,1) + do j = lbound(nest_var_grid,2), ubound(nest_var_grid,2) + var_grid((ioffset-1)*refine+i, (joffset-1)*refine+j) = nest_var_grid(i,j) + enddo + enddo + + end subroutine mn_overwrite_with_nest_init_values_r4 + + subroutine deallocate_land_mask_grids(land_mask_grids) + type(mn_land_mask_grids), intent(inout) :: land_mask_grids + + if (allocated(land_mask_grids%orog_grid)) deallocate(land_mask_grids%orog_grid) + if (allocated(land_mask_grids%orog_std_grid)) deallocate(land_mask_grids%orog_std_grid) + if (allocated(land_mask_grids%ls_mask_grid)) deallocate(land_mask_grids%ls_mask_grid) + if (allocated(land_mask_grids%soil_type_grid)) deallocate(land_mask_grids%soil_type_grid) + if (allocated(land_mask_grids%land_frac_grid)) deallocate(land_mask_grids%land_frac_grid) + if (allocated(land_mask_grids%geolat_grid)) deallocate(land_mask_grids%geolat_grid) + if (allocated(land_mask_grids%geolon_grid)) deallocate(land_mask_grids%geolon_grid) + end subroutine deallocate_land_mask_grids + + + + subroutine alloc_set_facwf(fix_grids) + type(mn_fix_grids), intent(inout) :: fix_grids + + integer :: i,j + + allocate(fix_grids%facwf_grid(lbound(fix_grids%facsf_grid,1):ubound(fix_grids%facsf_grid,1),lbound(fix_grids%facsf_grid,2):ubound(fix_grids%facsf_grid,2))) + + ! For land points, set facwf = 1.0 - facsf + ! To match initialization behavior, set any -999s to 0 + do i=lbound(fix_grids%facsf_grid,1),ubound(fix_grids%facsf_grid,1) + do j=lbound(fix_grids%facsf_grid,2),ubound(fix_grids%facsf_grid,2) + if (fix_grids%facsf_grid(i,j) .lt. -100) then + fix_grids%facsf_grid(i,j) = 0 + fix_grids%facwf_grid(i,j) = 0 + else + fix_grids%facwf_grid(i,j) = 1.0 - fix_grids%facsf_grid(i,j) + endif + enddo + enddo + end subroutine alloc_set_facwf + + + + subroutine mn_static_overwrite_ls_from_nest(fp_ls, nest_ls, refine, ioffset, joffset) + type(mn_land_mask_grids), intent(inout) :: fp_ls + type(mn_land_mask_grids), intent(in) :: nest_ls + integer, intent(in) :: refine, ioffset, joffset + + ! Update full panel with nest init values (there are a few mismatches) + ! TODO maybe add orog_raw/orog_filt + call mn_overwrite_with_nest_init_values("ls_mask", fp_ls%ls_mask_grid, nest_ls%ls_mask_grid, refine, ioffset, joffset) + + !if (is_fine_pe) then + ! call validate_navigation_fields("INIT", Atm_block, IPD_control, IPD_data, parent_grid_num, child_grid_num) + !endif + + call mn_overwrite_with_nest_init_values("soil_type", fp_ls%soil_type_grid, nest_ls%soil_type_grid, refine, ioffset, joffset) + call mn_overwrite_with_nest_init_values("land_frac", fp_ls%land_frac_grid, nest_ls%land_frac_grid, refine, ioffset, joffset) + + end subroutine mn_static_overwrite_ls_from_nest + + subroutine mn_static_overwrite_fix_from_nest(fp_fix, nest_fix, refine, ioffset, joffset) + type(mn_fix_grids), intent(inout) :: fp_fix + type(mn_fix_grids), intent(in) :: nest_fix + integer, intent(in) :: refine, ioffset, joffset + + call mn_overwrite_with_nest_init_values("deep_soil_temp", fp_fix%deep_soil_temp_grid, nest_fix%deep_soil_temp_grid, refine, ioffset, joffset) + call mn_overwrite_with_nest_init_values("veg_type", fp_fix%veg_type_grid, nest_fix%veg_type_grid, refine, ioffset, joffset) + call mn_overwrite_with_nest_init_values("slope_type", fp_fix%slope_type_grid, nest_fix%slope_type_grid, refine, ioffset, joffset) + call mn_overwrite_with_nest_init_values("max_snow_alb", fp_fix%max_snow_alb_grid, nest_fix%max_snow_alb_grid, refine, ioffset, joffset) + call mn_overwrite_with_nest_init_values("facsf", fp_fix%facsf_grid, nest_fix%facsf_grid, refine, ioffset, joffset) + call mn_overwrite_with_nest_init_values("facwf", fp_fix%facwf_grid, nest_fix%facwf_grid, refine, ioffset, joffset) + + call mn_overwrite_with_nest_init_values("alvsf", fp_fix%alvsf_grid, nest_fix%alvsf_grid, refine, ioffset, joffset) + call mn_overwrite_with_nest_init_values("alvwf", fp_fix%alvwf_grid, nest_fix%alvwf_grid, refine, ioffset, joffset) + call mn_overwrite_with_nest_init_values("alnsf", fp_fix%alnsf_grid, nest_fix%alnsf_grid, refine, ioffset, joffset) + call mn_overwrite_with_nest_init_values("alnwf", fp_fix%alnwf_grid, nest_fix%alnwf_grid, refine, ioffset, joffset) + + end subroutine mn_static_overwrite_fix_from_nest + + subroutine deallocate_fix_grids(fix_grids) + type(mn_fix_grids), intent(inout) :: fix_grids + + if (allocated(fix_grids%deep_soil_temp_grid)) deallocate(fix_grids%deep_soil_temp_grid) + if (allocated(fix_grids%veg_frac_grid)) deallocate(fix_grids%veg_frac_grid) + if (allocated(fix_grids%veg_type_grid)) deallocate(fix_grids%veg_type_grid) + !if (allocated(fix_grids%veg_greenness_grid)) deallocate(fix_grids%veg_greenness_grid) + if (allocated(fix_grids%slope_type_grid)) deallocate(fix_grids%slope_type_grid) + if (allocated(fix_grids%max_snow_alb_grid)) deallocate(fix_grids%max_snow_alb_grid) + if (allocated(fix_grids%facsf_grid)) deallocate(fix_grids%facsf_grid) + if (allocated(fix_grids%facwf_grid)) deallocate(fix_grids%facwf_grid) + if (allocated(fix_grids%alvsf_grid)) deallocate(fix_grids%alvsf_grid) + if (allocated(fix_grids%alvwf_grid)) deallocate(fix_grids%alvwf_grid) + if (allocated(fix_grids%alnsf_grid)) deallocate(fix_grids%alnsf_grid) + if (allocated(fix_grids%alnwf_grid)) deallocate(fix_grids%alnwf_grid) + + end subroutine deallocate_fix_grids + + + + + + + subroutine allocate_fv_moving_nest_prog_type(isd, ied, jsd, jed, npz, mn_prog) integer, intent(in) :: isd, ied, jsd, jed, npz type(fv_moving_nest_prog_type), intent(inout) :: mn_prog @@ -348,15 +537,17 @@ subroutine deallocate_fv_moving_nest_prog_type(mn_prog) end subroutine deallocate_fv_moving_nest_prog_type - subroutine allocate_fv_moving_nest_physics_type(isd, ied, jsd, jed, npz, move_physics, move_nsst, lsoil, nmtvr, levs, ntot2d, ntot3d, mn_phys) + subroutine allocate_fv_moving_nest_physics_type(isd, ied, jsd, jed, npz, move_physics, move_nsst, lsnow_lbound, lsnow_ubound, lsoil, nmtvr, levs, ntot2d, ntot3d, mn_phys) integer, intent(in) :: isd, ied, jsd, jed, npz logical, intent(in) :: move_physics, move_nsst - integer, intent(in) :: lsoil, nmtvr, levs, ntot2d, ntot3d ! From IPD_Control + integer, intent(in) :: lsnow_lbound, lsnow_ubound, lsoil, nmtvr, levs, ntot2d, ntot3d ! From IPD_Control type(fv_moving_nest_physics_type), intent(inout) :: mn_phys ! The local/temporary variables need to be allocated to the larger data (compute + halos) domain so that the nest motion code has halos to use allocate ( mn_phys%ts(isd:ied, jsd:jed) ) + !print '("[INFO] WDR allocate_fv_moving_nest_physics_type npe=",I0," lsnow_lbound=",I0," lsnow_ubound=",I0," lsoil=",I0)', mpp_pe(), lsnow_lbound, lsnow_ubound, lsoil + if (move_physics) then allocate ( mn_phys%slmsk(isd:ied, jsd:jed) ) allocate ( mn_phys%smc(isd:ied, jsd:jed, lsoil) ) @@ -441,6 +632,7 @@ subroutine allocate_fv_moving_nest_physics_type(isd, ied, jsd, jed, npz, move_p end if mn_phys%ts = +99999.9 + if (move_physics) then mn_phys%slmsk = +99999.9 mn_phys%smc = +99999.9 diff --git a/moving_nest/fv_moving_nest_utils.F90 b/moving_nest/fv_moving_nest_utils.F90 index 06136a2b1..d2795e4ca 100644 --- a/moving_nest/fv_moving_nest_utils.F90 +++ b/moving_nest/fv_moving_nest_utils.F90 @@ -66,13 +66,11 @@ module fv_moving_nest_utils_mod #else use IPD_typedefs, only: kind_phys => IPD_kind_phys #endif - #ifdef OVERLOAD_R4 use constantsR4_mod, only: grav #else use constants_mod, only: grav #endif - use boundary_mod, only: update_coarse_grid, update_coarse_grid_mpp use bounding_box_mod, only: bbox, bbox_get_C2F_index, fill_bbox use fms2_io_mod, only: read_data, write_data, open_file, close_file, register_axis, register_field @@ -94,6 +92,12 @@ module fv_moving_nest_utils_mod integer, parameter:: f_p = selected_real_kind(20) #endif +#ifdef OVERLOAD_R4 + real, parameter:: real_snan=x'FFBFFFFF' +#else + real, parameter:: real_snan=x'FFF7FFFFFFFFFFFF' +#endif + integer, parameter :: UWIND = 1 integer, parameter :: VWIND = 2 @@ -112,41 +116,50 @@ module fv_moving_nest_utils_mod interface fill_nest_halos_from_parent module procedure fill_nest_halos_from_parent_r4_2d - module procedure fill_nest_halos_from_parent_r4_3d + module procedure fill_nest_halos_from_parent_r4_3d_highz + module procedure fill_nest_halos_from_parent_r4_3d_lowhighz module procedure fill_nest_halos_from_parent_r4_4d module procedure fill_nest_halos_from_parent_r8_2d - module procedure fill_nest_halos_from_parent_r8_3d + module procedure fill_nest_halos_from_parent_r8_3d_highz + module procedure fill_nest_halos_from_parent_r8_3d_lowhighz module procedure fill_nest_halos_from_parent_r8_4d end interface fill_nest_halos_from_parent + interface alloc_halo_buffer module procedure alloc_halo_buffer_r4_2d - module procedure alloc_halo_buffer_r4_3d + module procedure alloc_halo_buffer_r4_3d_highz + module procedure alloc_halo_buffer_r4_3d_lowhighz module procedure alloc_halo_buffer_r4_4d module procedure alloc_halo_buffer_r8_2d - module procedure alloc_halo_buffer_r8_3d + module procedure alloc_halo_buffer_r8_3d_highz + module procedure alloc_halo_buffer_r8_3d_lowhighz module procedure alloc_halo_buffer_r8_4d end interface alloc_halo_buffer interface fill_nest_from_buffer module procedure fill_nest_from_buffer_r4_2d - module procedure fill_nest_from_buffer_r4_3d + module procedure fill_nest_from_buffer_r4_3d_highz + module procedure fill_nest_from_buffer_r4_3d_lowhighz module procedure fill_nest_from_buffer_r4_4d module procedure fill_nest_from_buffer_r8_2d - module procedure fill_nest_from_buffer_r8_3d + module procedure fill_nest_from_buffer_r8_3d_highz + module procedure fill_nest_from_buffer_r8_3d_lowhighz module procedure fill_nest_from_buffer_r8_4d end interface fill_nest_from_buffer interface fill_nest_from_buffer_cell_center module procedure fill_nest_from_buffer_cell_center_r4_2d - module procedure fill_nest_from_buffer_cell_center_r4_3d + module procedure fill_nest_from_buffer_cell_center_r4_3d_highz + module procedure fill_nest_from_buffer_cell_center_r4_3d_lowhighz module procedure fill_nest_from_buffer_cell_center_r4_4d module procedure fill_nest_from_buffer_cell_center_r8_2d - module procedure fill_nest_from_buffer_cell_center_r8_3d + module procedure fill_nest_from_buffer_cell_center_r8_3d_highz + module procedure fill_nest_from_buffer_cell_center_r8_3d_lowhighz module procedure fill_nest_from_buffer_cell_center_r8_4d end interface fill_nest_from_buffer_cell_center @@ -161,6 +174,28 @@ module fv_moving_nest_utils_mod module procedure fill_grid_from_supergrid_r8_4d end interface fill_grid_from_supergrid + ! Masked subroutines + interface fill_nest_halos_from_parent_masked + module procedure fill_nest_halos_from_parent_masked_r8_2d_const + module procedure fill_nest_halos_from_parent_masked_r8_2d_2d + module procedure fill_nest_halos_from_parent_masked_r8_3d_lowhighZ_const + module procedure fill_nest_halos_from_parent_masked_r8_3d_lowhighZ_1d + module procedure fill_nest_halos_from_parent_masked_r8_3d_lowhighZ_2d + end interface fill_nest_halos_from_parent_masked + + interface fill_nest_from_buffer_masked + module procedure fill_nest_from_buffer_masked_r8_2d_const + module procedure fill_nest_from_buffer_masked_r8_2d_2d + module procedure fill_nest_from_buffer_masked_r8_3d_1d + module procedure fill_nest_from_buffer_masked_r8_3d_2d + end interface fill_nest_from_buffer_masked + + interface fill_nest_from_buffer_cell_center_masked + module procedure fill_nest_from_buffer_cell_center_masked_2d_const + module procedure fill_nest_from_buffer_cell_center_masked_2d_2d + module procedure fill_nest_from_buffer_cell_center_masked_3d_1d + module procedure fill_nest_from_buffer_cell_center_masked_3d_2d + end interface fill_nest_from_buffer_cell_center_masked contains @@ -483,7 +518,7 @@ subroutine fill_nest_halos_from_parent_r8_2d(var_name, data_var, interp_type, wt end subroutine fill_nest_halos_from_parent_r8_2d - subroutine fill_nest_halos_from_parent_masked(var_name, data_var, interp_type, wt, ind, x_refine, y_refine, is_fine_pe, nest_domain, position, mask_var, mask_val, default_val) + subroutine fill_nest_halos_from_parent_masked_r8_2d_const(var_name, data_var, interp_type, wt, ind, x_refine, y_refine, is_fine_pe, nest_domain, position, mask_var, parent_mask_var, mask_val, default_val) character(len=*), intent(in) :: var_name real*8, allocatable, intent(inout) :: data_var(:,:) integer, intent(in) :: interp_type @@ -493,7 +528,8 @@ subroutine fill_nest_halos_from_parent_masked(var_name, data_var, interp_type, w logical, intent(in) :: is_fine_pe type(nest_domain_type), intent(inout) :: nest_domain integer, intent(in) :: position - real*4, allocatable, intent(in) :: mask_var(:,:) + real, allocatable, intent(in) :: mask_var(:,:) + real, allocatable, intent(in) :: parent_mask_var(:,:) integer, intent(in) :: mask_val real*8, intent(in) :: default_val @@ -529,10 +565,222 @@ subroutine fill_nest_halos_from_parent_masked(var_name, data_var, interp_type, w !! !!=========================================================== - call fill_nest_from_buffer_masked(interp_type, data_var, nbuffer, north_fine, north_coarse, NORTH, x_refine, y_refine, wt, ind, mask_var, mask_val, default_val) - call fill_nest_from_buffer_masked(interp_type, data_var, sbuffer, south_fine, south_coarse, SOUTH, x_refine, y_refine, wt, ind, mask_var, mask_val, default_val) - call fill_nest_from_buffer_masked(interp_type, data_var, ebuffer, east_fine, east_coarse, EAST, x_refine, y_refine, wt, ind, mask_var, mask_val, default_val) - call fill_nest_from_buffer_masked(interp_type, data_var, wbuffer, west_fine, west_coarse, WEST, x_refine, y_refine, wt, ind, mask_var, mask_val, default_val) + call fill_nest_from_buffer_masked(var_name, interp_type, data_var, nbuffer, north_fine, north_coarse, NORTH, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_val) + call fill_nest_from_buffer_masked(var_name, interp_type, data_var, sbuffer, south_fine, south_coarse, SOUTH, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_val) + call fill_nest_from_buffer_masked(var_name, interp_type, data_var, ebuffer, east_fine, east_coarse, EAST, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_val) + call fill_nest_from_buffer_masked(var_name, interp_type, data_var, wbuffer, west_fine, west_coarse, WEST, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_val) + + endif + + deallocate(nbuffer) + deallocate(sbuffer) + deallocate(ebuffer) + deallocate(wbuffer) + + end subroutine fill_nest_halos_from_parent_masked_r8_2d_const + + subroutine fill_nest_halos_from_parent_masked_r8_2d_2d(var_name, data_var, interp_type, wt, ind, x_refine, y_refine, is_fine_pe, nest_domain, position, mask_var, parent_mask_var, mask_val, default_grid) + character(len=*), intent(in) :: var_name + real*8, allocatable, intent(inout) :: data_var(:,:) + integer, intent(in) :: interp_type + real, allocatable, intent(in) :: wt(:,:,:) ! TODO should this also be real*8? + integer, allocatable, intent(in) :: ind(:,:,:) + integer, intent(in) :: x_refine, y_refine + logical, intent(in) :: is_fine_pe + type(nest_domain_type), intent(inout) :: nest_domain + integer, intent(in) :: position + real, allocatable, intent(in) :: mask_var(:,:) + real, allocatable, intent(in) :: parent_mask_var(:,:) + integer, intent(in) :: mask_val + real, allocatable, intent(in) :: default_grid(:,:) + + real*8, dimension(:,:), allocatable :: nbuffer, sbuffer, ebuffer, wbuffer + type(bbox) :: north_fine, north_coarse + type(bbox) :: south_fine, south_coarse + type(bbox) :: east_fine, east_coarse + type(bbox) :: west_fine, west_coarse + integer :: this_pe + integer :: nest_level = 1 ! TODO allow to vary + + this_pe = mpp_pe() + + !!=========================================================== + !! + !! Fill halo buffers + !! + !!=========================================================== + + call alloc_halo_buffer(nbuffer, north_fine, north_coarse, nest_domain, NORTH, position) + call alloc_halo_buffer(sbuffer, south_fine, south_coarse, nest_domain, SOUTH, position) + call alloc_halo_buffer(ebuffer, east_fine, east_coarse, nest_domain, EAST, position) + call alloc_halo_buffer(wbuffer, west_fine, west_coarse, nest_domain, WEST, position) + + ! Passes data from coarse grid to fine grid's halo + call mpp_update_nest_fine(data_var, nest_domain, wbuffer, sbuffer, ebuffer, nbuffer, nest_level, position=position) + + if (is_fine_pe) then + + !!=========================================================== + !! + !! Apply halo data + !! + !!=========================================================== + + call fill_nest_from_buffer_masked(var_name, interp_type, data_var, nbuffer, north_fine, north_coarse, NORTH, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_grid) + call fill_nest_from_buffer_masked(var_name, interp_type, data_var, sbuffer, south_fine, south_coarse, SOUTH, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_grid) + call fill_nest_from_buffer_masked(var_name, interp_type, data_var, ebuffer, east_fine, east_coarse, EAST, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_grid) + call fill_nest_from_buffer_masked(var_name, interp_type, data_var, wbuffer, west_fine, west_coarse, WEST, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_grid) + + endif + + deallocate(nbuffer) + deallocate(sbuffer) + deallocate(ebuffer) + deallocate(wbuffer) + + end subroutine fill_nest_halos_from_parent_masked_r8_2d_2d + + + + + subroutine fill_nest_halos_from_parent_masked_r8_3d_lowhighz_const(var_name, data_var, interp_type, wt, ind, x_refine, y_refine, is_fine_pe, nest_domain, position, low_z, high_z, mask_var, parent_mask_var, mask_val, default_val) + character(len=*), intent(in) :: var_name + real*8, allocatable, intent(inout) :: data_var(:,:,:) + integer, intent(in) :: interp_type + real, allocatable, intent(in) :: wt(:,:,:) ! TODO should this be real*8? + integer, allocatable, intent(in) :: ind(:,:,:) + integer, intent(in) :: x_refine, y_refine + logical, intent(in) :: is_fine_pe + type(nest_domain_type), intent(inout) :: nest_domain + integer, intent(in) :: position, low_z, high_z + real, allocatable, intent(in) :: mask_var(:,:) + real, allocatable, intent(in) :: parent_mask_var(:,:) + integer, intent(in) :: mask_val + real*8, intent(in) :: default_val + + real*8 :: default_vector(low_z:high_z) + + default_vector = default_val + + call fill_nest_halos_from_parent_masked_r8_3d_lowhighz_1d(var_name, data_var, interp_type, wt, ind, x_refine, y_refine, is_fine_pe, nest_domain, position, low_z, high_z, mask_var, parent_mask_var, mask_val, default_vector) + + end subroutine fill_nest_halos_from_parent_masked_r8_3d_lowhighz_const + + + + subroutine fill_nest_halos_from_parent_masked_r8_3d_lowhighz_1d(var_name, data_var, interp_type, wt, ind, x_refine, y_refine, is_fine_pe, nest_domain, position, low_z, high_z, mask_var, parent_mask_var, mask_val, default_val) + character(len=*), intent(in) :: var_name + real*8, allocatable, intent(inout) :: data_var(:,:,:) + integer, intent(in) :: interp_type + real, allocatable, intent(in) :: wt(:,:,:) ! TODO should this be real*8? + integer, allocatable, intent(in) :: ind(:,:,:) + integer, intent(in) :: x_refine, y_refine + logical, intent(in) :: is_fine_pe + type(nest_domain_type), intent(inout) :: nest_domain + integer, intent(in) :: position, low_z, high_z + real, allocatable, intent(in) :: mask_var(:,:) + real, allocatable, intent(in) :: parent_mask_var(:,:) + integer, intent(in) :: mask_val + real*8, intent(in) :: default_val(low_z:high_z) + + real*8, dimension(:,:,:), allocatable :: nbuffer, sbuffer, ebuffer, wbuffer + type(bbox) :: north_fine, north_coarse + type(bbox) :: south_fine, south_coarse + type(bbox) :: east_fine, east_coarse + type(bbox) :: west_fine, west_coarse + integer :: this_pe + integer :: nest_level = 1 ! TODO allow to vary + + this_pe = mpp_pe() + + !!=========================================================== + !! + !! Fill halo buffers + !! + !!=========================================================== + + call alloc_halo_buffer(nbuffer, north_fine, north_coarse, nest_domain, NORTH, position, low_z, high_z) + call alloc_halo_buffer(sbuffer, south_fine, south_coarse, nest_domain, SOUTH, position, low_z, high_z) + call alloc_halo_buffer(ebuffer, east_fine, east_coarse, nest_domain, EAST, position, low_z, high_z) + call alloc_halo_buffer(wbuffer, west_fine, west_coarse, nest_domain, WEST, position, low_z, high_z) + + ! Passes data from coarse grid to fine grid's halo + call mpp_update_nest_fine(data_var, nest_domain, wbuffer, sbuffer, ebuffer, nbuffer, nest_level, position=position) + + if (is_fine_pe) then + + !!=========================================================== + !! + !! Apply halo data + !! + !!=========================================================== + + call fill_nest_from_buffer_masked(var_name, interp_type, data_var, nbuffer, north_fine, north_coarse, low_z, high_z, NORTH, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_val) + call fill_nest_from_buffer_masked(var_name, interp_type, data_var, sbuffer, south_fine, south_coarse, low_z, high_z, SOUTH, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_val) + call fill_nest_from_buffer_masked(var_name, interp_type, data_var, ebuffer, east_fine, east_coarse, low_z, high_z, EAST, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_val) + call fill_nest_from_buffer_masked(var_name, interp_type, data_var, wbuffer, west_fine, west_coarse, low_z, high_z, WEST, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_val) + + endif + + deallocate(nbuffer) + deallocate(sbuffer) + deallocate(ebuffer) + deallocate(wbuffer) + + end subroutine fill_nest_halos_from_parent_masked_r8_3d_lowhighz_1d + + + subroutine fill_nest_halos_from_parent_masked_r8_3d_lowhighz_2d(var_name, data_var, interp_type, wt, ind, x_refine, y_refine, is_fine_pe, nest_domain, position, low_z, high_z, mask_var, parent_mask_var, mask_val, default_grid) + character(len=*), intent(in) :: var_name + real*8, allocatable, intent(inout) :: data_var(:,:,:) + integer, intent(in) :: interp_type + real, allocatable, intent(in) :: wt(:,:,:) ! TODO should this be real*8? + integer, allocatable, intent(in) :: ind(:,:,:) + integer, intent(in) :: x_refine, y_refine + logical, intent(in) :: is_fine_pe + type(nest_domain_type), intent(inout) :: nest_domain + integer, intent(in) :: position, low_z, high_z + real, allocatable, intent(in) :: mask_var(:,:) + real, allocatable, intent(in) :: parent_mask_var(:,:) + integer, intent(in) :: mask_val + real, allocatable, intent(in) :: default_grid(:,:) + + real*8, dimension(:,:,:), allocatable :: nbuffer, sbuffer, ebuffer, wbuffer + type(bbox) :: north_fine, north_coarse + type(bbox) :: south_fine, south_coarse + type(bbox) :: east_fine, east_coarse + type(bbox) :: west_fine, west_coarse + integer :: this_pe + integer :: nest_level = 1 ! TODO allow to vary + + this_pe = mpp_pe() + + !!=========================================================== + !! + !! Fill halo buffers + !! + !!=========================================================== + + call alloc_halo_buffer(nbuffer, north_fine, north_coarse, nest_domain, NORTH, position, low_z, high_z) + call alloc_halo_buffer(sbuffer, south_fine, south_coarse, nest_domain, SOUTH, position, low_z, high_z) + call alloc_halo_buffer(ebuffer, east_fine, east_coarse, nest_domain, EAST, position, low_z, high_z) + call alloc_halo_buffer(wbuffer, west_fine, west_coarse, nest_domain, WEST, position, low_z, high_z) + + ! Passes data from coarse grid to fine grid's halo + call mpp_update_nest_fine(data_var, nest_domain, wbuffer, sbuffer, ebuffer, nbuffer, nest_level, position=position) + + if (is_fine_pe) then + + !!=========================================================== + !! + !! Apply halo data + !! + !!=========================================================== + + call fill_nest_from_buffer_masked(var_name, interp_type, data_var, nbuffer, north_fine, north_coarse, low_z, high_z, NORTH, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_grid) + call fill_nest_from_buffer_masked(var_name, interp_type, data_var, sbuffer, south_fine, south_coarse, low_z, high_z, SOUTH, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_grid) + call fill_nest_from_buffer_masked(var_name, interp_type, data_var, ebuffer, east_fine, east_coarse, low_z, high_z, EAST, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_grid) + call fill_nest_from_buffer_masked(var_name, interp_type, data_var, wbuffer, west_fine, west_coarse, low_z, high_z, WEST, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_grid) endif @@ -541,10 +789,11 @@ subroutine fill_nest_halos_from_parent_masked(var_name, data_var, interp_type, w deallocate(ebuffer) deallocate(wbuffer) - end subroutine fill_nest_halos_from_parent_masked + end subroutine fill_nest_halos_from_parent_masked_r8_3d_lowhighz_2d + - subroutine fill_nest_halos_from_parent_r4_3d(var_name, data_var, interp_type, wt, ind, x_refine, y_refine, is_fine_pe, nest_domain, position, nz) + subroutine fill_nest_halos_from_parent_r4_3d_highz(var_name, data_var, interp_type, wt, ind, x_refine, y_refine, is_fine_pe, nest_domain, position, nz) character(len=*), intent(in) :: var_name real*4, allocatable, intent(inout) :: data_var(:,:,:) integer, intent(in) :: interp_type @@ -555,6 +804,22 @@ subroutine fill_nest_halos_from_parent_r4_3d(var_name, data_var, interp_type, wt type(nest_domain_type), intent(inout) :: nest_domain integer, intent(in) :: position, nz + + call fill_nest_halos_from_parent_r4_3d_lowhighz(var_name, data_var, interp_type, wt, ind, x_refine, y_refine, is_fine_pe, nest_domain, position, 1, nz) + + end subroutine fill_nest_halos_from_parent_r4_3d_highz + + subroutine fill_nest_halos_from_parent_r4_3d_lowhighz(var_name, data_var, interp_type, wt, ind, x_refine, y_refine, is_fine_pe, nest_domain, position, low_z, high_z) + character(len=*), intent(in) :: var_name + real*4, allocatable, intent(inout) :: data_var(:,:,:) + integer, intent(in) :: interp_type + real, allocatable, intent(in) :: wt(:,:,:) + integer, allocatable, intent(in) :: ind(:,:,:) + integer, intent(in) :: x_refine, y_refine + logical, intent(in) :: is_fine_pe + type(nest_domain_type), intent(inout) :: nest_domain + integer, intent(in) :: position, low_z, high_z + real*4, dimension(:,:,:), allocatable :: nbuffer, sbuffer, ebuffer, wbuffer type(bbox) :: north_fine, north_coarse type(bbox) :: south_fine, south_coarse @@ -571,10 +836,10 @@ subroutine fill_nest_halos_from_parent_r4_3d(var_name, data_var, interp_type, wt !! !!=========================================================== - call alloc_halo_buffer(nbuffer, north_fine, north_coarse, nest_domain, NORTH, position, nz) - call alloc_halo_buffer(sbuffer, south_fine, south_coarse, nest_domain, SOUTH, position, nz) - call alloc_halo_buffer(ebuffer, east_fine, east_coarse, nest_domain, EAST, position, nz) - call alloc_halo_buffer(wbuffer, west_fine, west_coarse, nest_domain, WEST, position, nz) + call alloc_halo_buffer(nbuffer, north_fine, north_coarse, nest_domain, NORTH, position, low_z, high_z) + call alloc_halo_buffer(sbuffer, south_fine, south_coarse, nest_domain, SOUTH, position, low_z, high_z) + call alloc_halo_buffer(ebuffer, east_fine, east_coarse, nest_domain, EAST, position, low_z, high_z) + call alloc_halo_buffer(wbuffer, west_fine, west_coarse, nest_domain, WEST, position, low_z, high_z) ! Passes data from coarse grid to fine grid's halo call mpp_update_nest_fine(data_var, nest_domain, wbuffer, sbuffer, ebuffer, nbuffer, nest_level, position=position) @@ -587,10 +852,10 @@ subroutine fill_nest_halos_from_parent_r4_3d(var_name, data_var, interp_type, wt !! !!=========================================================== - call fill_nest_from_buffer(interp_type, data_var, nbuffer, north_fine, north_coarse, nz, NORTH, x_refine, y_refine, wt, ind) - call fill_nest_from_buffer(interp_type, data_var, sbuffer, south_fine, south_coarse, nz, SOUTH, x_refine, y_refine, wt, ind) - call fill_nest_from_buffer(interp_type, data_var, ebuffer, east_fine, east_coarse, nz, EAST, x_refine, y_refine, wt, ind) - call fill_nest_from_buffer(interp_type, data_var, wbuffer, west_fine, west_coarse, nz, WEST, x_refine, y_refine, wt, ind) + call fill_nest_from_buffer(interp_type, data_var, nbuffer, north_fine, north_coarse, low_z, high_z, NORTH, x_refine, y_refine, wt, ind) + call fill_nest_from_buffer(interp_type, data_var, sbuffer, south_fine, south_coarse, low_z, high_z, SOUTH, x_refine, y_refine, wt, ind) + call fill_nest_from_buffer(interp_type, data_var, ebuffer, east_fine, east_coarse, low_z, high_z, EAST, x_refine, y_refine, wt, ind) + call fill_nest_from_buffer(interp_type, data_var, wbuffer, west_fine, west_coarse, low_z, high_z, WEST, x_refine, y_refine, wt, ind) endif @@ -599,10 +864,9 @@ subroutine fill_nest_halos_from_parent_r4_3d(var_name, data_var, interp_type, wt deallocate(ebuffer) deallocate(wbuffer) - end subroutine fill_nest_halos_from_parent_r4_3d + end subroutine fill_nest_halos_from_parent_r4_3d_lowhighz - - subroutine fill_nest_halos_from_parent_r8_3d(var_name, data_var, interp_type, wt, ind, x_refine, y_refine, is_fine_pe, nest_domain, position, nz) + subroutine fill_nest_halos_from_parent_r8_3d_highz(var_name, data_var, interp_type, wt, ind, x_refine, y_refine, is_fine_pe, nest_domain, position, nz) character(len=*), intent(in) :: var_name real*8, allocatable, intent(inout) :: data_var(:,:,:) integer, intent(in) :: interp_type @@ -613,6 +877,21 @@ subroutine fill_nest_halos_from_parent_r8_3d(var_name, data_var, interp_type, wt type(nest_domain_type), intent(inout) :: nest_domain integer, intent(in) :: position, nz + call fill_nest_halos_from_parent_r8_3d_lowhighz(var_name, data_var, interp_type, wt, ind, x_refine, y_refine, is_fine_pe, nest_domain, position, 1, nz) + + end subroutine fill_nest_halos_from_parent_r8_3d_highz + + subroutine fill_nest_halos_from_parent_r8_3d_lowhighz(var_name, data_var, interp_type, wt, ind, x_refine, y_refine, is_fine_pe, nest_domain, position, low_z, high_z) + character(len=*), intent(in) :: var_name + real*8, allocatable, intent(inout) :: data_var(:,:,:) + integer, intent(in) :: interp_type + real, allocatable, intent(in) :: wt(:,:,:) ! TODO should this be real*8? + integer, allocatable, intent(in) :: ind(:,:,:) + integer, intent(in) :: x_refine, y_refine + logical, intent(in) :: is_fine_pe + type(nest_domain_type), intent(inout) :: nest_domain + integer, intent(in) :: position, low_z, high_z + real*8, dimension(:,:,:), allocatable :: nbuffer, sbuffer, ebuffer, wbuffer type(bbox) :: north_fine, north_coarse type(bbox) :: south_fine, south_coarse @@ -629,10 +908,10 @@ subroutine fill_nest_halos_from_parent_r8_3d(var_name, data_var, interp_type, wt !! !!=========================================================== - call alloc_halo_buffer(nbuffer, north_fine, north_coarse, nest_domain, NORTH, position, nz) - call alloc_halo_buffer(sbuffer, south_fine, south_coarse, nest_domain, SOUTH, position, nz) - call alloc_halo_buffer(ebuffer, east_fine, east_coarse, nest_domain, EAST, position, nz) - call alloc_halo_buffer(wbuffer, west_fine, west_coarse, nest_domain, WEST, position, nz) + call alloc_halo_buffer(nbuffer, north_fine, north_coarse, nest_domain, NORTH, position, low_z, high_z) + call alloc_halo_buffer(sbuffer, south_fine, south_coarse, nest_domain, SOUTH, position, low_z, high_z) + call alloc_halo_buffer(ebuffer, east_fine, east_coarse, nest_domain, EAST, position, low_z, high_z) + call alloc_halo_buffer(wbuffer, west_fine, west_coarse, nest_domain, WEST, position, low_z, high_z) ! Passes data from coarse grid to fine grid's halo call mpp_update_nest_fine(data_var, nest_domain, wbuffer, sbuffer, ebuffer, nbuffer, nest_level, position=position) @@ -645,10 +924,10 @@ subroutine fill_nest_halos_from_parent_r8_3d(var_name, data_var, interp_type, wt !! !!=========================================================== - call fill_nest_from_buffer(interp_type, data_var, nbuffer, north_fine, north_coarse, nz, NORTH, x_refine, y_refine, wt, ind) - call fill_nest_from_buffer(interp_type, data_var, sbuffer, south_fine, south_coarse, nz, SOUTH, x_refine, y_refine, wt, ind) - call fill_nest_from_buffer(interp_type, data_var, ebuffer, east_fine, east_coarse, nz, EAST, x_refine, y_refine, wt, ind) - call fill_nest_from_buffer(interp_type, data_var, wbuffer, west_fine, west_coarse, nz, WEST, x_refine, y_refine, wt, ind) + call fill_nest_from_buffer(interp_type, data_var, nbuffer, north_fine, north_coarse, low_z, high_z, NORTH, x_refine, y_refine, wt, ind) + call fill_nest_from_buffer(interp_type, data_var, sbuffer, south_fine, south_coarse, low_z, high_z, SOUTH, x_refine, y_refine, wt, ind) + call fill_nest_from_buffer(interp_type, data_var, ebuffer, east_fine, east_coarse, low_z, high_z, EAST, x_refine, y_refine, wt, ind) + call fill_nest_from_buffer(interp_type, data_var, wbuffer, west_fine, west_coarse, low_z, high_z, WEST, x_refine, y_refine, wt, ind) endif @@ -657,7 +936,7 @@ subroutine fill_nest_halos_from_parent_r8_3d(var_name, data_var, interp_type, wt deallocate(ebuffer) deallocate(wbuffer) - end subroutine fill_nest_halos_from_parent_r8_3d + end subroutine fill_nest_halos_from_parent_r8_3d_lowhighz subroutine fill_nest_halos_from_parent_r4_4d(var_name, data_var, interp_type, wt, ind, x_refine, y_refine, is_fine_pe, nest_domain, position, nz) @@ -836,17 +1115,27 @@ subroutine alloc_halo_buffer_r4_2d(buffer, bbox_fine, bbox_coarse, nest_domain, end subroutine alloc_halo_buffer_r4_2d - subroutine alloc_halo_buffer_r4_3d(buffer, bbox_fine, bbox_coarse, nest_domain, direction, position, nz) + subroutine alloc_halo_buffer_r4_3d_highz(buffer, bbox_fine, bbox_coarse, nest_domain, direction, position, high_z) real*4, dimension(:,:,:), allocatable, intent(out) :: buffer type(bbox), intent(out) :: bbox_fine, bbox_coarse type(nest_domain_type), intent(in) :: nest_domain - integer, intent(in) :: direction, position, nz + integer, intent(in) :: direction, position, high_z + + call alloc_halo_buffer_r4_3d_lowhighz(buffer, bbox_fine, bbox_coarse, nest_domain, direction, position, 1, high_z) + + end subroutine alloc_halo_buffer_r4_3d_highz + + subroutine alloc_halo_buffer_r4_3d_lowhighz(buffer, bbox_fine, bbox_coarse, nest_domain, direction, position, low_z, high_z) + real*4, dimension(:,:,:), allocatable, intent(out) :: buffer + type(bbox), intent(out) :: bbox_fine, bbox_coarse + type(nest_domain_type), intent(in) :: nest_domain + integer, intent(in) :: direction, position, low_z, high_z call bbox_get_C2F_index(nest_domain, bbox_fine, bbox_coarse, direction, position) if( bbox_coarse.ie .GE. bbox_coarse.is .AND. bbox_coarse.je .GE. bbox_coarse.js ) then - allocate(buffer(bbox_coarse.is:bbox_coarse.ie, bbox_coarse.js:bbox_coarse.je,1:nz)) + allocate(buffer(bbox_coarse.is:bbox_coarse.ie, bbox_coarse.js:bbox_coarse.je, low_z:high_z)) else ! The buffer must have some storage allocated, whether it's a useful buffer or just a dummy. allocate(buffer(1,1,1)) @@ -854,19 +1143,29 @@ subroutine alloc_halo_buffer_r4_3d(buffer, bbox_fine, bbox_coarse, nest_domain, buffer = 0 - end subroutine alloc_halo_buffer_r4_3d + end subroutine alloc_halo_buffer_r4_3d_lowhighz + + + subroutine alloc_halo_buffer_r8_3d_highz(buffer, bbox_fine, bbox_coarse, nest_domain, direction, position, high_z) + real*8, dimension(:,:,:), allocatable, intent(out) :: buffer + type(bbox), intent(out) :: bbox_fine, bbox_coarse + type(nest_domain_type), intent(in) :: nest_domain + integer, intent(in) :: direction, position, high_z + + call alloc_halo_buffer_r8_3d_lowhighz(buffer, bbox_fine, bbox_coarse, nest_domain, direction, position, 1, high_z) + end subroutine alloc_halo_buffer_r8_3d_highz - subroutine alloc_halo_buffer_r8_3d(buffer, bbox_fine, bbox_coarse, nest_domain, direction, position, nz) + subroutine alloc_halo_buffer_r8_3d_lowhighz(buffer, bbox_fine, bbox_coarse, nest_domain, direction, position, low_z, high_z) real*8, dimension(:,:,:), allocatable, intent(out) :: buffer type(bbox), intent(out) :: bbox_fine, bbox_coarse type(nest_domain_type), intent(in) :: nest_domain - integer, intent(in) :: direction, position, nz + integer, intent(in) :: direction, position, low_z, high_z call bbox_get_C2F_index(nest_domain, bbox_fine, bbox_coarse, direction, position) if( bbox_coarse.ie .GE. bbox_coarse.is .AND. bbox_coarse.je .GE. bbox_coarse.js ) then - allocate(buffer(bbox_coarse.is:bbox_coarse.ie, bbox_coarse.js:bbox_coarse.je,1:nz)) + allocate(buffer(bbox_coarse.is:bbox_coarse.ie, bbox_coarse.js:bbox_coarse.je, low_z:high_z)) else ! The buffer must have some storage allocated, whether it's a useful buffer or just a dummy. allocate(buffer(1,1,1)) @@ -874,7 +1173,7 @@ subroutine alloc_halo_buffer_r8_3d(buffer, bbox_fine, bbox_coarse, nest_domain, buffer = 0 - end subroutine alloc_halo_buffer_r8_3d + end subroutine alloc_halo_buffer_r8_3d_lowhighz subroutine alloc_halo_buffer_r4_4d(buffer, bbox_fine, bbox_coarse, nest_domain, direction, position, nz, n4d) @@ -1522,9 +1821,10 @@ subroutine fill_nest_from_buffer_r8_2d(interp_type, x, buffer, bbox_fine, bbox_c end subroutine fill_nest_from_buffer_r8_2d - subroutine fill_nest_from_buffer_masked(interp_type, x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind, mask_var, mask_val, default_val) + subroutine fill_nest_from_buffer_masked_r8_2d_const(var_name, interp_type, x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_val) implicit none + character(len=*), intent(in) :: var_name integer, intent(in) :: interp_type real*8, allocatable, intent(inout) :: x(:,:) real*8, allocatable, intent(in) :: buffer(:,:) @@ -1533,6 +1833,7 @@ subroutine fill_nest_from_buffer_masked(interp_type, x, buffer, bbox_fine, bbox_ real, allocatable, intent(in) :: wt(:,:,:) ! The final dimension is always 4 integer, allocatable, intent(in) :: ind(:,:,:) real, allocatable, intent(in) :: mask_var(:,:) + real, allocatable, intent(in) :: parent_mask_var(:,:) integer, intent(in) :: mask_val real*8, intent(in) :: default_val @@ -1542,12 +1843,14 @@ subroutine fill_nest_from_buffer_masked(interp_type, x, buffer, bbox_fine, bbox_ ! Output the interpolation type select case (interp_type) case (1) + print '("[WARN] fv_moving_nest_utils.F90 fill_nest_from_buffer_mask interp_type 1 not implemented. var_name=",A16)', var_name call fill_nest_from_buffer_cell_center("A", x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind) ! case (3) ! C grid staggered case (4) + print '("[WARN] fv_moving_nest_utils.F90 fill_nest_from_buffer_mask interp_type 4 not implemented. var_name=",A16)', var_name call fill_nest_from_buffer_cell_center("D", x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind) case (7) - call fill_nest_from_buffer_cell_center_masked("A", x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind, mask_var, mask_val, default_val) + call fill_nest_from_buffer_cell_center_masked(var_name, "A", x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_val) case (9) !call fill_nest_from_buffer_nearest_neighbor(x, buffer, bbox_fine, bbox_coarse, dir, wt) call mpp_error(FATAL, '2D fill_nest_from_buffer_nearest_neighbor not yet implemented.') @@ -1555,21 +1858,23 @@ subroutine fill_nest_from_buffer_masked(interp_type, x, buffer, bbox_fine, bbox_ call mpp_error(FATAL, 'interp_single_nest got invalid value for interp_type from namelist.') end select - end subroutine fill_nest_from_buffer_masked - + end subroutine fill_nest_from_buffer_masked_r8_2d_const - - subroutine fill_nest_from_buffer_r4_3d(interp_type, x, buffer, bbox_fine, bbox_coarse, nz, dir, x_refine, y_refine, wt, ind) + subroutine fill_nest_from_buffer_masked_r8_2d_2d(var_name, interp_type, x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_grid) implicit none + character(len=*), intent(in) :: var_name integer, intent(in) :: interp_type - real*4, allocatable, intent(inout) :: x(:,:,:) - real*4, allocatable, intent(in) :: buffer(:,:,:) + real*8, allocatable, intent(inout) :: x(:,:) + real*8, allocatable, intent(in) :: buffer(:,:) type(bbox), intent(in) :: bbox_fine, bbox_coarse - integer, intent(in) :: nz integer, intent(in) :: dir, x_refine, y_refine real, allocatable, intent(in) :: wt(:,:,:) ! The final dimension is always 4 integer, allocatable, intent(in) :: ind(:,:,:) + real, allocatable, intent(in) :: mask_var(:,:) + real, allocatable, intent(in) :: parent_mask_var(:,:) + integer, intent(in) :: mask_val + real, allocatable, intent(in) :: default_grid(:,:) integer :: this_pe this_pe = mpp_pe() @@ -1577,31 +1882,39 @@ subroutine fill_nest_from_buffer_r4_3d(interp_type, x, buffer, bbox_fine, bbox_c ! Output the interpolation type select case (interp_type) case (1) - call fill_nest_from_buffer_cell_center("A", x, buffer, bbox_fine, bbox_coarse, nz, dir, x_refine, y_refine, wt, ind) - ! case (3) ! C grid staggered + print '("[WARN] fv_moving_nest_utils.F90 fill_nest_from_buffer_mask interp_type 1 not implemented. var_name=",A16)', var_name + call fill_nest_from_buffer_cell_center("A", x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind) + ! case (3) ! C grid staggered case (4) - call fill_nest_from_buffer_cell_center("D", x, buffer, bbox_fine, bbox_coarse, nz, dir, x_refine, y_refine, wt, ind) + print '("[WARN] fv_moving_nest_utils.F90 fill_nest_from_buffer_mask interp_type 4 not implemented. var_name=",A16)', var_name + call fill_nest_from_buffer_cell_center("D", x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind) + case (7) + call fill_nest_from_buffer_cell_center_masked(var_name, "A", x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_grid) case (9) - !call fill_nest_from_buffer_nearest_neighbor(x, buffer, bbox_fine, bbox_coarse, nz, dir, wt) - call mpp_error(FATAL, 'fill_nest_from_buffer_nearest_neighbor is not yet implemented.') + !call fill_nest_from_buffer_nearest_neighbor(x, buffer, bbox_fine, bbox_coarse, dir, wt) + call mpp_error(FATAL, '2D fill_nest_from_buffer_nearest_neighbor not yet implemented.') case default call mpp_error(FATAL, 'interp_single_nest got invalid value for interp_type from namelist.') end select - end subroutine fill_nest_from_buffer_r4_3d + end subroutine fill_nest_from_buffer_masked_r8_2d_2d - subroutine fill_nest_from_buffer_r8_3d(interp_type, x, buffer, bbox_fine, bbox_coarse, nz, dir, x_refine, y_refine, wt, ind) + subroutine fill_nest_from_buffer_masked_r8_3d_1d(var_name, interp_type, x, buffer, bbox_fine, bbox_coarse, low_z, high_z, dir, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_vector) implicit none + character(len=*), intent(in) :: var_name integer, intent(in) :: interp_type real*8, allocatable, intent(inout) :: x(:,:,:) real*8, allocatable, intent(in) :: buffer(:,:,:) type(bbox), intent(in) :: bbox_fine, bbox_coarse - integer, intent(in) :: nz - integer, intent(in) :: dir, x_refine, y_refine + integer, intent(in) :: dir, x_refine, y_refine, low_z, high_z real, allocatable, intent(in) :: wt(:,:,:) ! The final dimension is always 4 integer, allocatable, intent(in) :: ind(:,:,:) + real, allocatable, intent(in) :: mask_var(:,:) + real, allocatable, intent(in) :: parent_mask_var(:,:) + integer, intent(in) :: mask_val + real*8, intent(in) :: default_vector(low_z:high_z) integer :: this_pe this_pe = mpp_pe() @@ -1609,22 +1922,167 @@ subroutine fill_nest_from_buffer_r8_3d(interp_type, x, buffer, bbox_fine, bbox_c ! Output the interpolation type select case (interp_type) case (1) - call fill_nest_from_buffer_cell_center("A", x, buffer, bbox_fine, bbox_coarse, nz, dir, x_refine, y_refine, wt, ind) + print '("[WARN] fv_moving_nest_utils.F90 fill_nest_from_buffer_mask interp_type 1 not implemented. var_name=",A16)', var_name + call mpp_error(FATAL, '3D fill_nest_from_buffer_nearest_neighbor not yet implemented.') + !call fill_nest_from_buffer_cell_center("A", x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind) ! case (3) ! C grid staggered case (4) - call fill_nest_from_buffer_cell_center("D", x, buffer, bbox_fine, bbox_coarse, nz, dir, x_refine, y_refine, wt, ind) + print '("[WARN] fv_moving_nest_utils.F90 fill_nest_from_buffer_mask interp_type 4 not implemented. var_name=",A16)', var_name + call mpp_error(FATAL, '3D fill_nest_from_buffer_nearest_neighbor not yet implemented.') + !call fill_nest_from_buffer_cell_center("D", x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind) + case (7) + call fill_nest_from_buffer_cell_center_masked(var_name, "A", x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, low_z, high_z, default_vector) case (9) - call mpp_error(FATAL, 'nearest_neighbor is not yet implemented for fv_moving_nest_utils.F90::fill_nest_from_buffer_3D_kindphys') - !call fill_nest_from_buffer_nearest_neighbor(x, buffer, bbox_fine, bbox_coarse, nz, dir, wt) + !call fill_nest_from_buffer_nearest_neighbor(x, buffer, bbox_fine, bbox_coarse, dir, wt) + call mpp_error(FATAL, '3D fill_nest_from_buffer_nearest_neighbor not yet implemented.') case default call mpp_error(FATAL, 'interp_single_nest got invalid value for interp_type from namelist.') end select - end subroutine fill_nest_from_buffer_r8_3d + end subroutine fill_nest_from_buffer_masked_r8_3d_1d - !>@brief This subroutine fills the nest halo data from the coarse grid data by downscaling. - !>@details Applicable to any interpolation type + subroutine fill_nest_from_buffer_masked_r8_3d_2d(var_name, interp_type, x, buffer, bbox_fine, bbox_coarse, low_z, high_z, dir, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_grid) + implicit none + + character(len=*), intent(in) :: var_name + integer, intent(in) :: interp_type + real*8, allocatable, intent(inout) :: x(:,:,:) + real*8, allocatable, intent(in) :: buffer(:,:,:) + type(bbox), intent(in) :: bbox_fine, bbox_coarse + integer, intent(in) :: dir, x_refine, y_refine, low_z, high_z + real, allocatable, intent(in) :: wt(:,:,:) ! The final dimension is always 4 + integer, allocatable, intent(in) :: ind(:,:,:) + real, allocatable, intent(in) :: mask_var(:,:) + real, allocatable, intent(in) :: parent_mask_var(:,:) + integer, intent(in) :: mask_val + real, allocatable, intent(in) :: default_grid(:,:) + + integer :: this_pe + this_pe = mpp_pe() + + ! Output the interpolation type + select case (interp_type) + case (1) + print '("[WARN] fv_moving_nest_utils.F90 fill_nest_from_buffer_mask interp_type 1 not implemented. var_name=",A16)', var_name + call mpp_error(FATAL, '3D fill_nest_from_buffer_nearest_neighbor not yet implemented.') + !call fill_nest_from_buffer_cell_center("A", x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind) + ! case (3) ! C grid staggered + case (4) + print '("[WARN] fv_moving_nest_utils.F90 fill_nest_from_buffer_mask interp_type 4 not implemented. var_name=",A16)', var_name + call mpp_error(FATAL, '3D fill_nest_from_buffer_nearest_neighbor not yet implemented.') + !call fill_nest_from_buffer_cell_center("D", x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind) + case (7) + call fill_nest_from_buffer_cell_center_masked(var_name, "A", x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, low_z, high_z, default_grid) + case (9) + !call fill_nest_from_buffer_nearest_neighbor(x, buffer, bbox_fine, bbox_coarse, dir, wt) + call mpp_error(FATAL, '3D fill_nest_from_buffer_nearest_neighbor not yet implemented.') + case default + call mpp_error(FATAL, 'interp_single_nest got invalid value for interp_type from namelist.') + end select + + end subroutine fill_nest_from_buffer_masked_r8_3d_2d + + + + subroutine fill_nest_from_buffer_r4_3d_highz(interp_type, x, buffer, bbox_fine, bbox_coarse, nz, dir, x_refine, y_refine, wt, ind) + implicit none + + integer, intent(in) :: interp_type + real*4, allocatable, intent(inout) :: x(:,:,:) + real*4, allocatable, intent(in) :: buffer(:,:,:) + type(bbox), intent(in) :: bbox_fine, bbox_coarse + integer, intent(in) :: nz + integer, intent(in) :: dir, x_refine, y_refine + real, allocatable, intent(in) :: wt(:,:,:) ! The final dimension is always 4 + integer, allocatable, intent(in) :: ind(:,:,:) + + call fill_nest_from_buffer_r4_3d_lowhighz(interp_type, x, buffer, bbox_fine, bbox_coarse, 1, nz, dir, x_refine, y_refine, wt, ind) + + end subroutine fill_nest_from_buffer_r4_3d_highz + + subroutine fill_nest_from_buffer_r4_3d_lowhighz(interp_type, x, buffer, bbox_fine, bbox_coarse, low_z, high_z, dir, x_refine, y_refine, wt, ind) + implicit none + + integer, intent(in) :: interp_type + real*4, allocatable, intent(inout) :: x(:,:,:) + real*4, allocatable, intent(in) :: buffer(:,:,:) + type(bbox), intent(in) :: bbox_fine, bbox_coarse + integer, intent(in) :: low_z, high_z + integer, intent(in) :: dir, x_refine, y_refine + real, allocatable, intent(in) :: wt(:,:,:) ! The final dimension is always 4 + integer, allocatable, intent(in) :: ind(:,:,:) + + integer :: this_pe + this_pe = mpp_pe() + + ! Output the interpolation type + select case (interp_type) + case (1) + call fill_nest_from_buffer_cell_center("A", x, buffer, bbox_fine, bbox_coarse, low_z, high_z, dir, x_refine, y_refine, wt, ind) + ! case (3) ! C grid staggered + case (4) + call fill_nest_from_buffer_cell_center("D", x, buffer, bbox_fine, bbox_coarse, low_z, high_z, dir, x_refine, y_refine, wt, ind) + case (9) + !call fill_nest_from_buffer_nearest_neighbor(x, buffer, bbox_fine, bbox_coarse, low_z, high_z, dir, wt) + call mpp_error(FATAL, 'fill_nest_from_buffer_nearest_neighbor is not yet implemented.') + case default + call mpp_error(FATAL, 'interp_single_nest got invalid value for interp_type from namelist.') + end select + + end subroutine fill_nest_from_buffer_r4_3d_lowhighz + + + subroutine fill_nest_from_buffer_r8_3d_highz(interp_type, x, buffer, bbox_fine, bbox_coarse, nz, dir, x_refine, y_refine, wt, ind) + implicit none + + integer, intent(in) :: interp_type + real*8, allocatable, intent(inout) :: x(:,:,:) + real*8, allocatable, intent(in) :: buffer(:,:,:) + type(bbox), intent(in) :: bbox_fine, bbox_coarse + integer, intent(in) :: nz + integer, intent(in) :: dir, x_refine, y_refine + real, allocatable, intent(in) :: wt(:,:,:) ! The final dimension is always 4 + integer, allocatable, intent(in) :: ind(:,:,:) + + call fill_nest_from_buffer_r8_3d_lowhighz(interp_type, x, buffer, bbox_fine, bbox_coarse, 1, nz, dir, x_refine, y_refine, wt, ind) + + end subroutine fill_nest_from_buffer_r8_3d_highz + + subroutine fill_nest_from_buffer_r8_3d_lowhighz(interp_type, x, buffer, bbox_fine, bbox_coarse, low_z, high_z, dir, x_refine, y_refine, wt, ind) + implicit none + + integer, intent(in) :: interp_type + real*8, allocatable, intent(inout) :: x(:,:,:) + real*8, allocatable, intent(in) :: buffer(:,:,:) + type(bbox), intent(in) :: bbox_fine, bbox_coarse + integer, intent(in) :: low_z, high_z + integer, intent(in) :: dir, x_refine, y_refine + real, allocatable, intent(in) :: wt(:,:,:) ! The final dimension is always 4 + integer, allocatable, intent(in) :: ind(:,:,:) + + integer :: this_pe + this_pe = mpp_pe() + + ! Output the interpolation type + select case (interp_type) + case (1) + call fill_nest_from_buffer_cell_center("A", x, buffer, bbox_fine, bbox_coarse, low_z, high_z, dir, x_refine, y_refine, wt, ind) + ! case (3) ! C grid staggered + case (4) + call fill_nest_from_buffer_cell_center("D", x, buffer, bbox_fine, bbox_coarse, low_z, high_z, dir, x_refine, y_refine, wt, ind) + case (9) + call mpp_error(FATAL, 'nearest_neighbor is not yet implemented for fv_moving_nest_utils.F90::fill_nest_from_buffer_3D_kindphys') + !call fill_nest_from_buffer_nearest_neighbor(x, buffer, bbox_fine, bbox_coarse, low_z, high_z, dir, wt) + case default + call mpp_error(FATAL, 'interp_single_nest got invalid value for interp_type from namelist.') + end select + + end subroutine fill_nest_from_buffer_r8_3d_lowhighz + + + !>@brief This subroutine fills the nest halo data from the coarse grid data by downscaling. + !>@details Applicable to any interpolation type subroutine fill_nest_from_buffer_r4_4d(interp_type, x, buffer, bbox_fine, bbox_coarse, nz, dir, x_refine, y_refine, wt, ind) implicit none @@ -1795,8 +2253,9 @@ subroutine fill_nest_from_buffer_cell_center_r8_2d(stagger, x, buffer, bbox_fine end subroutine fill_nest_from_buffer_cell_center_r8_2d - subroutine fill_nest_from_buffer_cell_center_masked(stagger, x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind, mask_var, mask_val, default_val) + subroutine fill_nest_from_buffer_cell_center_masked_2d_const(var_name, stagger, x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_val) implicit none + character(len=*), intent(in) :: var_name character ( len = 1 ), intent(in) :: stagger real*8, allocatable, intent(inout) :: x(:,:) real*8, allocatable, intent(in) :: buffer(:,:) @@ -1805,12 +2264,21 @@ subroutine fill_nest_from_buffer_cell_center_masked(stagger, x, buffer, bbox_fin real, allocatable, intent(in) :: wt(:,:,:) ! The final dimension is always 4 integer, allocatable, intent(in) :: ind(:,:,:) real, allocatable, intent(in) :: mask_var(:,:) + real, allocatable, intent(in) :: parent_mask_var(:,:) integer, intent(in) :: mask_val real*8, intent(in) :: default_val character(len=8) :: dir_str integer :: i, j, k, ic, jc real :: tw + real :: dummy_val, dummy_mask + integer :: num_reset, num_weights + integer :: this_pe + + this_pe = mpp_pe() + + num_reset = 0 + dummy_val = real_snan select case(dir) case (NORTH) @@ -1842,19 +2310,74 @@ subroutine fill_nest_from_buffer_cell_center_masked(stagger, x, buffer, bbox_fin !if (mask_var(i,j) .eq. mask_val) then x(i,j) = 0.0 tw = 0.0 - if (buffer(ic,jc) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,1)*buffer(ic, jc ) - if (buffer(ic,jc+1) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,1)*buffer(ic, jc+1) - if (buffer(ic+1,jc+1) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,1)*buffer(ic+1,jc+1) - if (buffer(ic+1,jc) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,1)*buffer(ic+1,jc ) + num_weights = 0 + +! WDR Original -- seems like the wt values should range from 1-4, not all use wt(i,j,1) +! will likely alter land values of shifted physics fields in regression tests. +! old values were (slightly) incorrect -- averaged of the 4 nearby points instead of actual weights +! if (buffer(ic,jc) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,1)*buffer(ic, jc ) +! if (buffer(ic,jc+1) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,1)*buffer(ic, jc+1) +! if (buffer(ic+1,jc+1) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,1)*buffer(ic+1,jc+1) +! if (buffer(ic+1,jc) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,1)*buffer(ic+1,jc ) +! +! if (buffer(ic,jc) .gt. -1.0) tw = tw + wt(i,j,1) +! if (buffer(ic,jc+1) .gt. -1.0) tw = tw + wt(i,j,1) +! if (buffer(ic+1,jc+1) .gt. -1.0) tw = tw + wt(i,j,1) +! if (buffer(ic+1,jc) .gt. -1.0) tw = tw + wt(i,j,1) + + +! Intermediate: Corrected the weights +! if (buffer(ic,jc) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,1)*buffer(ic, jc ) +! if (buffer(ic,jc+1) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,2)*buffer(ic, jc+1) +! if (buffer(ic+1,jc+1) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,3)*buffer(ic+1,jc+1) +! if (buffer(ic+1,jc) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,4)*buffer(ic+1,jc ) +! +! if (buffer(ic,jc) .gt. -1.0) tw = tw + wt(i,j,1) +! if (buffer(ic,jc+1) .gt. -1.0) tw = tw + wt(i,j,2) +! if (buffer(ic+1,jc+1) .gt. -1.0) tw = tw + wt(i,j,3) +! if (buffer(ic+1,jc) .gt. -1.0) tw = tw + wt(i,j,4) + + +! print '("[INFO] MASK2D npe=",I0," ",A16," parent_mask_var(",I0,",",I0,")=",F15.5," mask_var(",I0,",",I0,")=",F15.5)', mpp_pe(), var_name, ic, jc, parent_mask_var(ic,jc), i, j, mask_var(i,j) + + + ! Note that weights don't seem to always be exactly 0.0 when the corner points are aligned + ! Use the land sea mask to choose which points to add to weight and buffer + if (parent_mask_var(ic,jc) .eq. mask_var(i,j) .and. wt(i,j,1) .gt. 0.0001 ) then + num_weights = num_weights + 1 + x(i,j) = x(i,j) + wt(i,j,1)*buffer(ic, jc ) + tw = tw + wt(i,j,1) + + endif + + if (parent_mask_var(ic,jc+1) .eq. mask_var(i,j) .and. wt(i,j,2) .gt. 0.0001) then + num_weights = num_weights + 2 + x(i,j) = x(i,j) + wt(i,j,2)*buffer(ic, jc+1) + tw = tw + wt(i,j,2) + + endif + + if (parent_mask_var(ic+1,jc+1) .eq. mask_var(i,j) .and. wt(i,j,3) .gt. 0.0001) then + num_weights = num_weights + 4 + x(i,j) = x(i,j) + wt(i,j,3)*buffer(ic+1,jc+1) + tw = tw + wt(i,j,3) + + endif + + if (parent_mask_var(ic+1,jc) .eq. mask_var(i,j) .and. wt(i,j,4) .gt. 0.0001) then + num_weights = num_weights + 8 + x(i,j) = x(i,j) + wt(i,j,4)*buffer(ic+1,jc ) + tw = tw + wt(i,j,4) + + endif - if (buffer(ic,jc) .gt. -1.0) tw = tw + wt(i,j,1) - if (buffer(ic,jc+1) .gt. -1.0) tw = tw + wt(i,j,1) - if (buffer(ic+1,jc+1) .gt. -1.0) tw = tw + wt(i,j,1) - if (buffer(ic+1,jc) .gt. -1.0) tw = tw + wt(i,j,1) if (tw .gt. 0.0) then x(i,j) = x(i,j) / tw else + num_reset = num_reset + 1 + dummy_val = buffer(ic, jc) + dummy_mask = mask_var(i,j) x(i,j) = default_val endif @@ -1862,10 +2385,357 @@ subroutine fill_nest_from_buffer_cell_center_masked(stagger, x, buffer, bbox_fin enddo endif - end subroutine fill_nest_from_buffer_cell_center_masked +! if (.not. isnan(dummy_val)) print '("[INFO] WDR fill_nest_from_buffer_cell_center_masked npe=",I0," num_reset=",I0," var=",A12," mask_var=",F10.4," dummy_val=",F14.4," ",E15.8)', mpp_pe(), num_reset, trim(var_name), dummy_mask, dummy_val, dummy_val + + end subroutine fill_nest_from_buffer_cell_center_masked_2d_const + + subroutine fill_nest_from_buffer_cell_center_masked_2d_2d(var_name, stagger, x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, default_grid) + implicit none + character(len=*), intent(in) :: var_name + character ( len = 1 ), intent(in) :: stagger + real*8, allocatable, intent(inout) :: x(:,:) + real*8, allocatable, intent(in) :: buffer(:,:) + type(bbox), intent(in) :: bbox_fine, bbox_coarse + integer, intent(in) :: dir, x_refine, y_refine + real, allocatable, intent(in) :: wt(:,:,:) ! The final dimension is always 4 + integer, allocatable, intent(in) :: ind(:,:,:) + real, allocatable, intent(in) :: mask_var(:,:) + real, allocatable, intent(in) :: parent_mask_var(:,:) + integer, intent(in) :: mask_val + real, allocatable, intent(in) :: default_grid(:,:) + + character(len=8) :: dir_str + integer :: i, j, k, ic, jc + real :: tw + real :: dummy_val, dummy_mask + integer :: num_reset, num_weights + integer :: this_pe + + this_pe = mpp_pe() + + num_reset = 0 + dummy_val = real_snan + + select case(dir) + case (NORTH) + dir_str = "NORTH" + case (SOUTH) + dir_str = "SOUTH" + case (EAST) + dir_str = "EAST" + case (WEST) + dir_str = "WEST" + case default + dir_str = "ERR DIR" + end select + + if( bbox_coarse%ie .GE. bbox_coarse%is .AND. bbox_coarse%je .GE. bbox_coarse%js ) then + do j=bbox_fine%js, bbox_fine%je + do i=bbox_fine%is, bbox_fine%ie + + ic = ind(i,j,1) + jc = ind(i,j,2) + + !x(i,j) = & + ! wt(i,j,1)*buffer(ic, jc ) + & + ! wt(i,j,2)*buffer(ic, jc+1) + & + ! wt(i,j,3)*buffer(ic+1,jc+1) + & + ! wt(i,j,4)*buffer(ic+1,jc ) + + ! Land type + !if (mask_var(i,j) .eq. mask_val) then + x(i,j) = 0.0 + tw = 0.0 + num_weights = 0 + +! WDR Original -- seems like the wt values should range from 1-4, not all use wt(i,j,1) +! will likely alter land values of shifted physics fields in regression tests. +! old values were (slightly) incorrect -- averaged of the 4 nearby points instead of actual weights +! if (buffer(ic,jc) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,1)*buffer(ic, jc ) +! if (buffer(ic,jc+1) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,1)*buffer(ic, jc+1) +! if (buffer(ic+1,jc+1) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,1)*buffer(ic+1,jc+1) +! if (buffer(ic+1,jc) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,1)*buffer(ic+1,jc ) +! +! if (buffer(ic,jc) .gt. -1.0) tw = tw + wt(i,j,1) +! if (buffer(ic,jc+1) .gt. -1.0) tw = tw + wt(i,j,1) +! if (buffer(ic+1,jc+1) .gt. -1.0) tw = tw + wt(i,j,1) +! if (buffer(ic+1,jc) .gt. -1.0) tw = tw + wt(i,j,1) + + +! Intermediate: Corrected the weights +! if (buffer(ic,jc) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,1)*buffer(ic, jc ) +! if (buffer(ic,jc+1) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,2)*buffer(ic, jc+1) +! if (buffer(ic+1,jc+1) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,3)*buffer(ic+1,jc+1) +! if (buffer(ic+1,jc) .gt. -1.0) x(i,j) = x(i,j) + wt(i,j,4)*buffer(ic+1,jc ) +! +! if (buffer(ic,jc) .gt. -1.0) tw = tw + wt(i,j,1) +! if (buffer(ic,jc+1) .gt. -1.0) tw = tw + wt(i,j,2) +! if (buffer(ic+1,jc+1) .gt. -1.0) tw = tw + wt(i,j,3) +! if (buffer(ic+1,jc) .gt. -1.0) tw = tw + wt(i,j,4) + + +! print '("[INFO] MASK2D npe=",I0," ",A16," parent_mask_var(",I0,",",I0,")=",F15.5," mask_var(",I0,",",I0,")=",F15.5)', mpp_pe(), var_name, ic, jc, parent_mask_var(ic,jc), i, j, mask_var(i,j) + + + + ! Note that weights don't seem to always be exactly 0.0 when the corner points are aligned + ! Use the land sea mask to choose which points to add to weight and buffer + if (parent_mask_var(ic,jc) .eq. mask_var(i,j) .and. wt(i,j,1) .gt. 0.0001 ) then + num_weights = num_weights + 1 + x(i,j) = x(i,j) + wt(i,j,1)*buffer(ic, jc ) + tw = tw + wt(i,j,1) + + endif + + if (parent_mask_var(ic,jc+1) .eq. mask_var(i,j) .and. wt(i,j,2) .gt. 0.0001) then + num_weights = num_weights + 2 + x(i,j) = x(i,j) + wt(i,j,2)*buffer(ic, jc+1) + tw = tw + wt(i,j,2) + + endif + + if (parent_mask_var(ic+1,jc+1) .eq. mask_var(i,j) .and. wt(i,j,3) .gt. 0.0001) then + num_weights = num_weights + 4 + x(i,j) = x(i,j) + wt(i,j,3)*buffer(ic+1,jc+1) + tw = tw + wt(i,j,3) + + endif + + if (parent_mask_var(ic+1,jc) .eq. mask_var(i,j) .and. wt(i,j,4) .gt. 0.0001) then + num_weights = num_weights + 8 + x(i,j) = x(i,j) + wt(i,j,4)*buffer(ic+1,jc ) + tw = tw + wt(i,j,4) + + endif + + if (tw .gt. 0.0) then + x(i,j) = x(i,j) / tw + else + num_reset = num_reset + 1 + dummy_val = buffer(ic, jc) + dummy_mask = mask_var(i,j) + x(i,j) = default_grid(i,j) + endif + + enddo + enddo + endif + +! if (.not. isnan(dummy_val)) print '("[INFO] WDR fill_nest_from_buffer_cell_center_masked npe=",I0," num_reset=",I0," var=",A12," mask_var=",F10.4," dummy_val=",F14.4," ",E15.8)', mpp_pe(), num_reset, trim(var_name), dummy_mask, dummy_val, dummy_val + + end subroutine fill_nest_from_buffer_cell_center_masked_2d_2d - subroutine fill_nest_from_buffer_cell_center_r4_3d(stagger, x, buffer, bbox_fine, bbox_coarse, nz, dir, x_refine, y_refine, wt, ind) + subroutine fill_nest_from_buffer_cell_center_masked_3d_1d(var_name, stagger, x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, low_z, high_z, default_vector) + implicit none + character(len=*), intent(in) :: var_name + character ( len = 1 ), intent(in) :: stagger + real*8, allocatable, intent(inout) :: x(:,:,:) + real*8, allocatable, intent(in) :: buffer(:,:,:) + type(bbox), intent(in) :: bbox_fine, bbox_coarse + integer, intent(in) :: dir, x_refine, y_refine + real, allocatable, intent(in) :: wt(:,:,:) ! The final dimension is always 4 + integer, allocatable, intent(in) :: ind(:,:,:) + real, allocatable, intent(in) :: mask_var(:,:) + real, allocatable, intent(in) :: parent_mask_var(:,:) + integer, intent(in) :: mask_val, low_z, high_z + real*8, intent(in) :: default_vector(low_z:high_z) + + character(len=8) :: dir_str + integer :: i, j, k, ic, jc + real :: tw + real :: dummy_val, dummy_mask + + dummy_val = real_snan + + select case(dir) + case (NORTH) + dir_str = "NORTH" + case (SOUTH) + dir_str = "SOUTH" + case (EAST) + dir_str = "EAST" + case (WEST) + dir_str = "WEST" + case default + dir_str = "ERR DIR" + end select + + if( bbox_coarse%ie .GE. bbox_coarse%is .AND. bbox_coarse%je .GE. bbox_coarse%js ) then + do j=bbox_fine%js, bbox_fine%je + do i=bbox_fine%is, bbox_fine%ie + + ic = ind(i,j,1) + jc = ind(i,j,2) + + !x(i,j) = & + ! wt(i,j,1)*buffer(ic, jc ) + & + ! wt(i,j,2)*buffer(ic, jc+1) + & + ! wt(i,j,3)*buffer(ic+1,jc+1) + & + ! wt(i,j,4)*buffer(ic+1,jc ) + + ! Land type + !if (mask_var(i,j) .eq. mask_val) then + + + do k=lbound(x,3), ubound(x,3) + x(i,j,k) = 0.0 + tw = 0.0 + + + +! print '("[INFO] MASK3D npe=",I0," ",A16," parent_mask_var(",I0,",",I0,")=",F15.5," mask_var(",I0,",",I0,")=",F15.5)', mpp_pe(), var_name, ic, jc, parent_mask_var(ic,jc), i, j, mask_var(i,j) + + + ! Use the land sea mask to choose which points to add to weight and buffer + if (parent_mask_var(ic,jc) .eq. mask_var(i,j)) then + x(i,j,k) = x(i,j,k) + wt(i,j,1)*buffer(ic, jc ,k) + tw = tw + wt(i,j,1) + endif + + if (parent_mask_var(ic,jc+1) .eq. mask_var(i,j)) then + x(i,j,k) = x(i,j,k) + wt(i,j,2)*buffer(ic, jc+1,k) + tw = tw + wt(i,j,2) + endif + + if (parent_mask_var(ic+1,jc+1) .eq. mask_var(i,j)) then + x(i,j,k) = x(i,j,k) + wt(i,j,3)*buffer(ic+1,jc+1,k) + tw = tw + wt(i,j,3) + endif + + if (parent_mask_var(ic+1,jc) .eq. mask_var(i,j)) then + x(i,j,k) = x(i,j,k) + wt(i,j,4)*buffer(ic+1,jc ,k) + tw = tw + wt(i,j,4) + endif + + + if (tw .gt. 0.0) then + x(i,j,k) = x(i,j,k) / tw + else + dummy_val = buffer(ic, jc,k) + dummy_mask = mask_var(i,j) + x(i,j,k) = default_vector(k) + endif + + enddo + enddo + enddo + endif + +! if (.not. isnan(dummy_val)) then +! print '("[INFO WDR CCM3 fill_nest_from_buffer_cell_center_masked_3d npe=",I0)', mpp_pe() +! print '("[INFO] WDR CCM3 fill_nest_from_buffer_cell_center_masked npe=",I0," var=",A16," mask_var=",F10.4," dummy_val=",F14.4," ",E15.8)', mpp_pe(), trim(var_name), dummy_mask, dummy_val, dummy_val +! endif + + end subroutine fill_nest_from_buffer_cell_center_masked_3d_1d + + + subroutine fill_nest_from_buffer_cell_center_masked_3d_2d(var_name, stagger, x, buffer, bbox_fine, bbox_coarse, dir, x_refine, y_refine, wt, ind, mask_var, parent_mask_var, mask_val, low_z, high_z, default_grid) + implicit none + character(len=*), intent(in) :: var_name + character ( len = 1 ), intent(in) :: stagger + real*8, allocatable, intent(inout) :: x(:,:,:) + real*8, allocatable, intent(in) :: buffer(:,:,:) + type(bbox), intent(in) :: bbox_fine, bbox_coarse + integer, intent(in) :: dir, x_refine, y_refine + real, allocatable, intent(in) :: wt(:,:,:) ! The final dimension is always 4 + integer, allocatable, intent(in) :: ind(:,:,:) + real, allocatable, intent(in) :: mask_var(:,:) + real, allocatable, intent(in) :: parent_mask_var(:,:) + integer, intent(in) :: mask_val, low_z, high_z + real, allocatable, intent(in) :: default_grid(:,:) + + character(len=8) :: dir_str + integer :: i, j, k, ic, jc + real :: tw + real :: dummy_val, dummy_mask + + dummy_val = real_snan + + select case(dir) + case (NORTH) + dir_str = "NORTH" + case (SOUTH) + dir_str = "SOUTH" + case (EAST) + dir_str = "EAST" + case (WEST) + dir_str = "WEST" + case default + dir_str = "ERR DIR" + end select + + if( bbox_coarse%ie .GE. bbox_coarse%is .AND. bbox_coarse%je .GE. bbox_coarse%js ) then + do j=bbox_fine%js, bbox_fine%je + do i=bbox_fine%is, bbox_fine%ie + + ic = ind(i,j,1) + jc = ind(i,j,2) + + !x(i,j) = & + ! wt(i,j,1)*buffer(ic, jc ) + & + ! wt(i,j,2)*buffer(ic, jc+1) + & + ! wt(i,j,3)*buffer(ic+1,jc+1) + & + ! wt(i,j,4)*buffer(ic+1,jc ) + + ! Land type + !if (mask_var(i,j) .eq. mask_val) then + + + do k=lbound(x,3), ubound(x,3) + x(i,j,k) = 0.0 + tw = 0.0 + + + +! print '("[INFO] MASK3D npe=",I0," ",A16," parent_mask_var(",I0,",",I0,")=",F15.5," mask_var(",I0,",",I0,")=",F15.5)', mpp_pe(), var_name, ic, jc, parent_mask_var(ic,jc), i, j, mask_var(i,j) + + + ! Use the land sea mask to choose which points to add to weight and buffer + if (parent_mask_var(ic,jc) .eq. mask_var(i,j)) then + x(i,j,k) = x(i,j,k) + wt(i,j,1)*buffer(ic, jc ,k) + tw = tw + wt(i,j,1) + endif + + if (parent_mask_var(ic,jc+1) .eq. mask_var(i,j)) then + x(i,j,k) = x(i,j,k) + wt(i,j,2)*buffer(ic, jc+1,k) + tw = tw + wt(i,j,2) + endif + + if (parent_mask_var(ic+1,jc+1) .eq. mask_var(i,j)) then + x(i,j,k) = x(i,j,k) + wt(i,j,3)*buffer(ic+1,jc+1,k) + tw = tw + wt(i,j,3) + endif + + if (parent_mask_var(ic+1,jc) .eq. mask_var(i,j)) then + x(i,j,k) = x(i,j,k) + wt(i,j,4)*buffer(ic+1,jc ,k) + tw = tw + wt(i,j,4) + endif + + + if (tw .gt. 0.0) then + x(i,j,k) = x(i,j,k) / tw + else + dummy_val = buffer(ic, jc,k) + dummy_mask = mask_var(i,j) + x(i,j,k) = default_grid(i,j) + endif + + enddo + enddo + enddo + endif + +! if (.not. isnan(dummy_val)) then +! print '("[INFO WDR CCM3 fill_nest_from_buffer_cell_center_masked_3d npe=",I0)', mpp_pe() +! print '("[INFO] WDR CCM3 fill_nest_from_buffer_cell_center_masked npe=",I0," var=",A16," mask_var=",F10.4," dummy_val=",F14.4," ",E15.8)', mpp_pe(), trim(var_name), dummy_mask, dummy_val, dummy_val +! endif + + end subroutine fill_nest_from_buffer_cell_center_masked_3d_2d + + + + subroutine fill_nest_from_buffer_cell_center_r4_3d_highz(stagger, x, buffer, bbox_fine, bbox_coarse, nz, dir, x_refine, y_refine, wt, ind) implicit none character ( len = 1 ), intent(in) :: stagger real*4, allocatable, intent(inout) :: x(:,:,:) @@ -1876,6 +2746,21 @@ subroutine fill_nest_from_buffer_cell_center_r4_3d(stagger, x, buffer, bbox_fine real, allocatable, intent(in) :: wt(:,:,:) ! The final dimension is always 4 integer, allocatable, intent(in) :: ind(:,:,:) + call fill_nest_from_buffer_cell_center_r4_3d_lowhighz(stagger, x, buffer, bbox_fine, bbox_coarse, 1, nz, dir, x_refine, y_refine, wt, ind) + + end subroutine fill_nest_from_buffer_cell_center_r4_3d_highz + + subroutine fill_nest_from_buffer_cell_center_r4_3d_lowhighz(stagger, x, buffer, bbox_fine, bbox_coarse, low_z, high_z, dir, x_refine, y_refine, wt, ind) + implicit none + character ( len = 1 ), intent(in) :: stagger + real*4, allocatable, intent(inout) :: x(:,:,:) + real*4, allocatable, intent(in) :: buffer(:,:,:) + type(bbox), intent(in) :: bbox_fine, bbox_coarse + integer, intent(in) :: low_z, high_z + integer, intent(in) :: dir, x_refine, y_refine + real, allocatable, intent(in) :: wt(:,:,:) ! The final dimension is always 4 + integer, allocatable, intent(in) :: ind(:,:,:) + character(len=8) :: dir_str integer :: i, j, k, ic, jc @@ -1893,7 +2778,7 @@ subroutine fill_nest_from_buffer_cell_center_r4_3d(stagger, x, buffer, bbox_fine end select if( bbox_coarse%ie .GE. bbox_coarse%is .AND. bbox_coarse%je .GE. bbox_coarse%js ) then - do k=1,nz + do k=low_z, high_z do j=bbox_fine%js, bbox_fine%je do i=bbox_fine%is, bbox_fine%ie !if (stagger == "A") then @@ -1915,9 +2800,9 @@ subroutine fill_nest_from_buffer_cell_center_r4_3d(stagger, x, buffer, bbox_fine enddo endif - end subroutine fill_nest_from_buffer_cell_center_r4_3d + end subroutine fill_nest_from_buffer_cell_center_r4_3d_lowhighz - subroutine fill_nest_from_buffer_cell_center_r8_3d(stagger, x, buffer, bbox_fine, bbox_coarse, nz, dir, x_refine, y_refine, wt, ind) + subroutine fill_nest_from_buffer_cell_center_r8_3d_highz(stagger, x, buffer, bbox_fine, bbox_coarse, nz, dir, x_refine, y_refine, wt, ind) implicit none character ( len = 1 ), intent(in) :: stagger real*8, allocatable, intent(inout) :: x(:,:,:) @@ -1928,6 +2813,21 @@ subroutine fill_nest_from_buffer_cell_center_r8_3d(stagger, x, buffer, bbox_fine real, allocatable, intent(in) :: wt(:,:,:) ! The final dimension is always 4 integer, allocatable, intent(in) :: ind(:,:,:) + call fill_nest_from_buffer_cell_center_r8_3d_lowhighz(stagger, x, buffer, bbox_fine, bbox_coarse, 1, nz, dir, x_refine, y_refine, wt, ind) + + end subroutine fill_nest_from_buffer_cell_center_r8_3d_highz + + subroutine fill_nest_from_buffer_cell_center_r8_3d_lowhighz(stagger, x, buffer, bbox_fine, bbox_coarse, low_z, high_z, dir, x_refine, y_refine, wt, ind) + implicit none + character ( len = 1 ), intent(in) :: stagger + real*8, allocatable, intent(inout) :: x(:,:,:) + real*8, allocatable, intent(in) :: buffer(:,:,:) + type(bbox), intent(in) :: bbox_fine, bbox_coarse + integer, intent(in) :: low_z, high_z + integer, intent(in) :: dir, x_refine, y_refine + real, allocatable, intent(in) :: wt(:,:,:) ! The final dimension is always 4 + integer, allocatable, intent(in) :: ind(:,:,:) + character(len=8) :: dir_str integer :: i, j, k, ic, jc @@ -1945,7 +2845,7 @@ subroutine fill_nest_from_buffer_cell_center_r8_3d(stagger, x, buffer, bbox_fine end select if( bbox_coarse%ie .GE. bbox_coarse%is .AND. bbox_coarse%je .GE. bbox_coarse%js ) then - do k=1,nz + do k=low_z, high_z do j=bbox_fine%js, bbox_fine%je do i=bbox_fine%is, bbox_fine%ie !if (stagger == "A") then @@ -1966,7 +2866,7 @@ subroutine fill_nest_from_buffer_cell_center_r8_3d(stagger, x, buffer, bbox_fine enddo endif - end subroutine fill_nest_from_buffer_cell_center_r8_3d + end subroutine fill_nest_from_buffer_cell_center_r8_3d_lowhighz subroutine fill_nest_from_buffer_cell_center_r4_4d(stagger, x, buffer, bbox_fine, bbox_coarse, nz, dir, x_refine, y_refine, wt, ind)