From 885dc427ee0ff460fe59d3d32e300d58d48dbaaa Mon Sep 17 00:00:00 2001 From: Miles A Curry Date: Fri, 30 Oct 2020 16:33:21 +0000 Subject: [PATCH 1/6] Remove previous interpolation of snow albedo --- .../mpas_init_atm_static.F | 97 +------------------ 1 file changed, 2 insertions(+), 95 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_static.F b/src/core_init_atmosphere/mpas_init_atm_static.F index 021b36c96e..796a3ee9a4 100644 --- a/src/core_init_atmosphere/mpas_init_atm_static.F +++ b/src/core_init_atmosphere/mpas_init_atm_static.F @@ -829,108 +829,15 @@ subroutine init_atm_static(mesh, dims, configs) ! if (trim(config_maxsnowalbedo_data) == 'MODIS') then + geog_sub_path = 'maxsnowalb_modis/' + call mpas_log_write('Using MODIS 0.05-deg data for maximum snow albedo') if (supersample_fac > 1) then call mpas_log_write(' Dataset will be supersampled by a factor of $i', intArgs=(/supersample_fac/)) end if - nx = 1206 - ny = 1206 - nz = 1 - isigned = 1 - endian = 0 - wordsize = 2 - scalefactor = 0.01 - msgval = real(-999.0,kind=R4KIND)*real(0.01,kind=R4KIND) - fillval = 0.0 - allocate(rarray(nx,ny,nz)) - allocate(nhs(nCells)) - nhs(:) = 0 snoalb(:) = 0.0 - rarray_ptr = c_loc(rarray) - - start_lat = 90.0 - 0.05 * 0.5 / supersample_fac - start_lon = -180.0 + 0.05 * 0.5 / supersample_fac - geog_sub_path = 'maxsnowalb_modis/' - - do jTileStart = 1,02401,ny-6 - jTileEnd = jTileStart + ny - 1 - 6 - - do iTileStart=1,06001,nx-6 - iTileEnd = iTileStart + nx - 1 - 6 - write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(geog_data_path)//trim(geog_sub_path), & - iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd - call mpas_log_write(trim(fname)) - call mpas_f_to_c_string(fname, c_fname) - - call read_geogrid(c_fname,rarray_ptr,nx,ny,nz,isigned,endian, & - wordsize,istatus) - call init_atm_check_read_error(istatus, fname) - rarray(:,:,:) = rarray(:,:,:) * scalefactor - - iPoint = 1 - do j=supersample_fac * 3 + 1, supersample_fac * (ny-3) - do i=supersample_fac * 3 + 1, supersample_fac * (nx-3) - ii = (i - 1) / supersample_fac + 1 - jj = (j - 1) / supersample_fac + 1 - - lat_pt = start_lat - (supersample_fac*(jTileStart-1) + j - (supersample_fac*3+1)) * 0.05 / supersample_fac - lon_pt = start_lon + (supersample_fac*(iTileStart-1) + i - (supersample_fac*3+1)) * 0.05 / supersample_fac - lat_pt = lat_pt * PI / 180.0 - lon_pt = lon_pt * PI / 180.0 - - iPoint = nearest_cell(lat_pt,lon_pt,iPoint,nCells,maxEdges, & - nEdgesOnCell,cellsOnCell, & - latCell,lonCell) - if (rarray(ii,jj,1) /= msgval) then - - ! - ! This field only matters for land cells, and for all but the outermost boundary cells, - ! we can safely assume that the nearest model grid cell contains the pixel (else, a different - ! cell would be nearest) - ! - if (landmask(iPoint) == 1 .and. bdyMaskCell(iPoint) < nBdyLayers) then - snoalb(iPoint) = snoalb(iPoint) + rarray(ii,jj,1) - nhs(iPoint) = nhs(iPoint) + 1 - - ! For outermost land cells, additional work is needed to verify that the pixel - ! actually lies within the nearest cell - else if (landmask(iPoint) == 1) then - zPixel = sphere_radius * sin(lat_pt) ! Model cell coordinates assume a "full" sphere radius - xPixel = sphere_radius * cos(lon_pt) * cos(lat_pt) ! at this point, so we need to ues the same radius - yPixel = sphere_radius * sin(lon_pt) * cos(lat_pt) ! for source pixel coordinates - - if (mpas_in_cell(xPixel, yPixel, zPixel, xCell(iPoint), yCell(iPoint), zCell(iPoint), & - nEdgesOnCell(iPoint), verticesOnCell(:,iPoint), xVertex, yVertex, zVertex)) then - snoalb(iPoint) = snoalb(iPoint) + rarray(ii,jj,1) - nhs(iPoint) = nhs(iPoint) + 1 - end if - end if - end if - end do - end do - - end do - end do - - do iCell = 1,nCells - ! - ! Mismatches in land mask can lead to MPAS land points with no maximum snow albedo. - ! Ideally, we would perform a search for nearby valid albedos, but for now using - ! the fill value will at least allow the model to run. In general, the number of cells - ! to be treated in this way tends to be a very small fraction of the total number of cells. - ! - if (nhs(iCell) == 0) then - snoalb(iCell) = fillval - else - snoalb(iCell) = snoalb(iCell) / real(nhs(iCell)) - end if - snoalb(iCell) = 0.01_RKIND * snoalb(iCell) ! Convert from percent to fraction - end do - deallocate(rarray) - deallocate(nhs) - else if (trim(config_maxsnowalbedo_data) == 'NCEP') then call mpas_log_write('Using NCEP 1.0-deg data for maximum snow albedo') From 6860e63718722a71427f09b1774baa8bf56550d8 Mon Sep 17 00:00:00 2001 From: Miles A Curry Date: Fri, 30 Oct 2020 17:38:07 +0000 Subject: [PATCH 2/6] Update geotile manager to work with negative grid spacing Before this commit, the geotile manager assumed that the grid spacing values of a geotile dataset, dx, and dy, would always be positive; however, datasets may start at the northern most latitude or eastern most longitude and move southward, eastward or both. This commit updates the geotile manager to use negative dx and dy values. --- src/core_init_atmosphere/mpas_geotile_manager.F | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core_init_atmosphere/mpas_geotile_manager.F b/src/core_init_atmosphere/mpas_geotile_manager.F index 0aa3298153..fd8e4f12e2 100644 --- a/src/core_init_atmosphere/mpas_geotile_manager.F +++ b/src/core_init_atmosphere/mpas_geotile_manager.F @@ -274,8 +274,8 @@ function mpas_geotile_mgr_init(mgr, path) result(ierr) ! Calculate the total number of pixels in x dir ! NOTE: This calculation assumes that a dataset is a global dataset and may ! not work correctly for non-global datasets - mgr % pixel_nx = nint(360.0_RKIND / dx) - mgr % pixel_ny = nint(180.0_RKIND / dy) + mgr % pixel_nx = nint(360.0_RKIND / abs(dx)) + mgr % pixel_ny = nint(180.0_RKIND / abs(dy)) ! Calculate the number of tiles in the x, y directions ! NOTE: This calculation assumes that a dataset is a global dataset and may From de8af90815f25f82673b5937f918c34bc020a536 Mon Sep 17 00:00:00 2001 From: Miles A Curry Date: Fri, 30 Oct 2020 18:03:32 +0000 Subject: [PATCH 3/6] Initialize, finalize and push all modis snoalb onto the stack This commit initializes and finalizes the geotile manager for the modis max snow albedo interpolation and it also loads and pushes all tiles onto the stack and pushes them off. It does not do any interpolation. --- .../mpas_init_atm_static.F | 37 +++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/src/core_init_atmosphere/mpas_init_atm_static.F b/src/core_init_atmosphere/mpas_init_atm_static.F index 796a3ee9a4..512c7d4c43 100644 --- a/src/core_init_atmosphere/mpas_init_atm_static.F +++ b/src/core_init_atmosphere/mpas_init_atm_static.F @@ -836,8 +836,45 @@ subroutine init_atm_static(mesh, dims, configs) call mpas_log_write(' Dataset will be supersampled by a factor of $i', intArgs=(/supersample_fac/)) end if + ierr = mgr % init(trim(config_geog_data_path)//trim(geog_sub_path)) + if (ierr /= 0) then + call mpas_log_write('Error occurred when initializing the interpolation of snow albedo (snoalb)', & + messageType=MPAS_LOG_CRIT) + endif + snoalb(:) = 0.0 + do iCell = 1, nCells + + tile => null() + ierr = mgr % get_tile(latCell(iCell), lonCell(iCell), tile) + if (ierr /= 0 .or. .not. associated(tile)) then + call mpas_log_write('Could not get the tile that contained cell $i', intArgs=(/iCell/), messageType=MPAS_LOG_CRIT) + end if + + if (tile % is_processed) then + cycle + end if + + tile % is_processed = .True. + + ierr = mgr % push_tile(tile) + if (ierr /= 0) then + call mpas_log_write('Error pushing this tile onto the stack: '//trim(tile%fname), messageType=MPAS_LOG_CRIT) + end if + end do + + do while(.not. mgr % is_stack_empty()) + tile => mgr % pop_tile() + call mpas_log_write('Processing tile: '//trim(tile % fname)) + end do + + ierr = mgr % finalize() + if (ierr /= 0) then + call mpas_log_write('Error occurred when finalizing the interpolation of snow albedo (snoalb)', & + messageType=MPAS_LOG_CRIT) + endif + else if (trim(config_maxsnowalbedo_data) == 'NCEP') then call mpas_log_write('Using NCEP 1.0-deg data for maximum snow albedo') From 315a6a5d590c42744aface142d53495f4c2b8a0a Mon Sep 17 00:00:00 2001 From: Miles A Curry Date: Wed, 4 Nov 2020 22:32:09 +0000 Subject: [PATCH 4/6] Enable the geotile manager to supersample tiles This commit enables the supersampling of tiles by updating the geotile manager tile_to_latlon function to subdivide the pixels of a tile if a supersampling factor is given in the function call. The code that calls tile_to_latlon will need to pass in supersampled coordinates if supersample_fac > 1. This is currently needed for the interpolation of maximum snow albedo, but will also be needed for future supersampling when other datasets become coarser than grid spacing. --- .../mpas_geotile_manager.F | 30 ++++++++++++++----- 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/src/core_init_atmosphere/mpas_geotile_manager.F b/src/core_init_atmosphere/mpas_geotile_manager.F index fd8e4f12e2..4913c7bec7 100644 --- a/src/core_init_atmosphere/mpas_geotile_manager.F +++ b/src/core_init_atmosphere/mpas_geotile_manager.F @@ -682,16 +682,22 @@ end function mpas_geotile_mgr_gen_tile_name ! ! public subroutine mpas_geotile_mgr_tile_to_latlon => tile_to_latlon ! - !> \brief Translate a tile indices to latitude and longitude, + !> \brief Find the latitude, longitude location of a tile's pixels !> \author Miles A. Curry !> \date 02/20/2020 !> \details - !> Given a tile, translate the relative tile coordinates, i, j, of that - !> tile to a latitude and longitude coordinate. Upon success, lat, lon - !> will be in the range of -1/2 * pi to 1/2 * pi and 0 to 2.0 * pi, respectively. + !> Given a tile, translate the pixel coordinates, i, j to a corresponding latitude + !> and longitude coordinate. + !> + !> If supersample_fac is present each pixel will be subdivided into supersample_fac ^ 2 + !> sub-pixels. If supersample_fac is greater than 1, then the calling code will need + !> to pass in supersampled i and j coordinates. + !> + !> Upon success, lat, lon will be in the range of -1/2 * pi to 1/2 * pi and 0 to + !> 2.0 * pi, respectively. ! !----------------------------------------------------------------------- - subroutine mpas_geotile_mgr_tile_to_latlon(mgr, tile, j, i, lat, lon) + subroutine mpas_geotile_mgr_tile_to_latlon(mgr, tile, j, i, lat, lon, supersample_fac) implicit none @@ -701,24 +707,34 @@ subroutine mpas_geotile_mgr_tile_to_latlon(mgr, tile, j, i, lat, lon) integer, value :: i real (kind=RKIND), intent(out) :: lat real (kind=RKIND), intent(out) :: lon + integer, optional, intent(in) :: supersample_fac integer, pointer :: tile_bdr real (kind=RKIND), pointer :: known_lon real (kind=RKIND), pointer :: known_lat real (kind=RKIND), pointer :: dx real (kind=RKIND), pointer :: dy + integer :: supersample_lcl integer :: ierr ierr = 0 + if (present(supersample_fac)) then + supersample_lcl = supersample_fac + else + supersample_lcl = 1 + end if + call mpas_pool_get_config(mgr % pool, 'tile_bdr', tile_bdr) call mpas_pool_get_config(mgr % pool, 'known_lat', known_lat) call mpas_pool_get_config(mgr % pool, 'known_lon', known_lon) call mpas_pool_get_config(mgr % pool, 'dx', dx) call mpas_pool_get_config(mgr % pool, 'dy', dy) - lat = real((j - (tile_bdr + 1) + tile % y - 1), kind=RKIND) * dy + known_lat - lon = real((i - (tile_bdr + 1) + tile % x - 1), kind=RKIND) * dx + known_lon + lat = known_lat + real(j - (supersample_lcl * tile_bdr + 1) + (supersample_lcl * (tile % y - 1)), kind=RKIND) * dy & + / supersample_lcl + lon = known_lon + real(i - (supersample_lcl * tile_bdr + 1) + (supersample_lcl * (tile % x - 1)), kind=RKIND) * dx & + / supersample_lcl call deg2Rad(lat) call deg2Rad(lon) From a8ab6f9686f117ac0df8224da17dfa10ebab62ac Mon Sep 17 00:00:00 2001 From: Miles A Curry Date: Wed, 4 Nov 2020 22:48:52 +0000 Subject: [PATCH 5/6] Parallel interpolation of maximum snow albedo This commit adds logic to interpolate maximum snow albedo in parallel. Because other fields have not yet been updated, parallelization of maximum snow albedo is disabled to do the check within init_atm_cases. --- .../mpas_init_atm_static.F | 136 +++++++++++++++--- 1 file changed, 120 insertions(+), 16 deletions(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_static.F b/src/core_init_atmosphere/mpas_init_atm_static.F index 512c7d4c43..2fe3b28cbe 100644 --- a/src/core_init_atmosphere/mpas_init_atm_static.F +++ b/src/core_init_atmosphere/mpas_init_atm_static.F @@ -129,10 +129,12 @@ subroutine init_atm_static(mesh, dims, configs) integer (kind=I8KIND), dimension(:), pointer :: ter_integer real (kind=RKIND), dimension(:), pointer :: soiltemp real (kind=RKIND), dimension(:), pointer :: snoalb + integer (kind=I8KIND), dimension(:), pointer :: snoalb_integer real (kind=RKIND), dimension(:), pointer :: shdmin, shdmax real (kind=RKIND), dimension(:,:), pointer :: greenfrac real (kind=RKIND), dimension(:,:), pointer :: albedo12m real (kind=RKIND) :: msgval, fillval + real (kind=RKIND), pointer :: missing_value integer, pointer :: category_min, category_max integer, dimension(:), pointer :: lu_index integer, dimension(:), pointer :: soilcat_top @@ -150,6 +152,7 @@ subroutine init_atm_static(mesh, dims, configs) type (mpas_geotile_type), pointer :: tile real (kind=RKIND) :: tval + integer (kind=I8KIND) :: i8val integer, pointer :: tile_bdr integer, pointer :: tile_nx, tile_ny @@ -842,33 +845,134 @@ subroutine init_atm_static(mesh, dims, configs) messageType=MPAS_LOG_CRIT) endif + call mpas_pool_get_config(mgr % pool, 'tile_bdr', tile_bdr) + call mpas_pool_get_config(mgr % pool, 'tile_x', tile_nx) + call mpas_pool_get_config(mgr % pool, 'tile_y', tile_ny) + call mpas_pool_get_config(mgr % pool, 'missing_value', missing_value) + call mpas_pool_get_config(mgr % pool, 'scale_factor', scalefactor_ptr) + scalefactor = scalefactor_ptr + + allocate(nhs(nCells)) + allocate(snoalb_integer(nCells)) + snoalb_integer(:) = 0 snoalb(:) = 0.0 + nhs(:) = 0 + fillval = 0.0 do iCell = 1, nCells + if (nhs(iCell) == 0) then + tile => null() + ierr = mgr % get_tile(latCell(iCell), lonCell(iCell), tile) + if (ierr /= 0 .or. .not. associated(tile)) then + call mpas_log_write('Could not get tile that contained cell $i', intArgs=(/iCell/), messageType=MPAS_LOG_CRIT) + end if - tile => null() - ierr = mgr % get_tile(latCell(iCell), lonCell(iCell), tile) - if (ierr /= 0 .or. .not. associated(tile)) then - call mpas_log_write('Could not get the tile that contained cell $i', intArgs=(/iCell/), messageType=MPAS_LOG_CRIT) - end if + ierr = mgr % push_tile(tile) + if (ierr /= 0) then + call mpas_log_write("Error pushing this tile onto the stack: "//trim(tile%fname), messageType=MPAS_LOG_CRIT) + end if + end if - if (tile % is_processed) then - cycle - end if + do while (.not. mgr % is_stack_empty()) + tile => mgr % pop_tile() - tile % is_processed = .True. + if (tile % is_processed) then + cycle + end if - ierr = mgr % push_tile(tile) - if (ierr /= 0) then - call mpas_log_write('Error pushing this tile onto the stack: '//trim(tile%fname), messageType=MPAS_LOG_CRIT) - end if + call mpas_log_write('Processing tile: '//trim(tile % fname)) + + all_pixels_mapped_to_halo_cells = .true. + + do j = supersample_fac * tile_bdr + 1, supersample_fac * (tile_ny + tile_bdr), 1 + do i = supersample_fac * tile_bdr + 1, supersample_fac * (tile_nx + tile_bdr), 1 + + ii = (i - 1) / supersample_fac + 1 + jj = (j - 1) / supersample_fac + 1 + + i8val = int(tile % tile(ii, jj, 1), kind=I8KIND) + + call mgr % tile_to_latlon(tile, j, i, lat_pt, lon_pt, supersample_fac) + call mpas_latlon_to_xyz(xPixel, yPixel, zPixel, sphere_radius, lat_pt, lon_pt) + call mpas_kd_search(tree, (/xPixel, yPixel, zPixel/), res) + + if (bdyMaskCell(res % id) < nBdyLayers) then + ! + ! This field only matters for land cells, and for all but the outermost boundary cells, + ! we can safely assume that the nearest model grid cell contains the pixel (else, a different + ! cell would be nearest). + ! + ! Since values in i8val are not yet scaled, we can compare them to missing_value, which + ! also is not scaled, without scaling either value + if (landmask(res % id) == 1 .and. i8val /= int(missing_value, kind=I8KIND)) then + snoalb_integer(res % id) = snoalb_integer(res % id) + i8val + nhs(res % id) = nhs(res % id) + 1 + end if + + ! + ! When a pixel maps to a non-land cell or is a missing value, the values are not accumulated + ! above; however, these pixels may still reside in an owned cell, in which case we will still need + ! to push the tile's neighbors onto the stack for processing. + ! + if (res % id <= nCellsSolve) then + all_pixels_mapped_to_halo_cells = .false. + end if + ! For outermost cells, additional work is needed to verify that the pixel + ! actually lies within the nearest cell + else + if (mpas_in_cell(xPixel, yPixel, zPixel, xCell(res % id), yCell(res % id), zCell(res % id), & + nEdgesOnCell(res % id), verticesOnCell(:,res % id), xVertex, yVertex, zVertex)) then + + ! Since values in i8val are not yet scaled, we can compare them to missing_value, which + ! also is not scaled, without scaling either value + if (landmask(res % id) == 1 .and. i8val /= int(missing_value, kind=I8KIND)) then + snoalb_integer(res % id) = snoalb_integer(res % id) + i8val + nhs(res % id) = nhs(res % id) + 1 + end if + + ! + ! When a pixel maps to a non-land cell or is a missing value, the values are not accumulated + ! above; however, these pixels may still reside in an owned cell, in which case we will still need + ! to push the tile's neighbors onto the stack for processing. + ! + if (res % id <= nCellsSolve) then + all_pixels_mapped_to_halo_cells = .false. + end if + end if + end if + end do + end do + + tile % is_processed = .true. + + if (.not. all_pixels_mapped_to_halo_cells) then + ierr = mgr % push_neighbors(tile) + if (ierr /= 0) then + call mpas_log_write("Error pushing the tile neighbors of: "//trim(tile%fname), messageType=MPAS_LOG_CRIT) + end if + end if + end do end do - do while(.not. mgr % is_stack_empty()) - tile => mgr % pop_tile() - call mpas_log_write('Processing tile: '//trim(tile % fname)) + do iCell = 1, nCells + ! + ! Mismatches in land mask can lead to MPAS land points with no maximum snow albedo. + ! Ideally, we would perform a search for nearby valid albedos, but for now using + ! the fill value will at least allow the model to run. In general, the number of cells + ! to be treated in this way tends to be a very small fraction of the total number of cells. + ! + if (nhs(iCell) == 0) then + snoalb(iCell) = fillval + else + snoalb(iCell) = real(snoalb_integer(iCell), kind=R8KIND) / real(nhs(iCell), kind=R8KIND) + snoalb(iCell) = snoalb(iCell) * scalefactor + snoalb(iCell) = 0.01_RKIND * snoalb(iCell) ! Convert from percent to fraction + endif end do + deallocate(nhs) + deallocate(snoalb_integer) + ierr = mgr % finalize() if (ierr /= 0) then call mpas_log_write('Error occurred when finalizing the interpolation of snow albedo (snoalb)', & From 9f4dfbf32b55dcf906213fa291fb7e207bd577ef Mon Sep 17 00:00:00 2001 From: Miles A Curry Date: Mon, 30 Nov 2020 23:33:34 +0000 Subject: [PATCH 6/6] Change scalefactor from c_float to RKIND Because read_geogrid is no longer takes scalefactor as an argument, it is no longer necessary for it to be a c_float. --- src/core_init_atmosphere/mpas_init_atm_static.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core_init_atmosphere/mpas_init_atm_static.F b/src/core_init_atmosphere/mpas_init_atm_static.F index 2fe3b28cbe..41007badd8 100644 --- a/src/core_init_atmosphere/mpas_init_atm_static.F +++ b/src/core_init_atmosphere/mpas_init_atm_static.F @@ -88,7 +88,7 @@ subroutine init_atm_static(mesh, dims, configs) integer,dimension(:,:),allocatable:: ncat real(kind=RKIND), pointer :: scalefactor_ptr - real(kind=c_float):: scalefactor + real(kind=RKIND) :: scalefactor real(kind=c_float),dimension(:,:,:),pointer,contiguous :: rarray type(c_ptr) :: rarray_ptr