Skip to content

Commit

Permalink
Merge branch 'init_atmosphere/snow_albd_interp' into develop (PR #750)
Browse files Browse the repository at this point in the history 
This merge allows the init_atmosphere model to safely interpolate maximum snow
albedo in parallel with any mesh partition file (not just CVT partition files).
However, because other fields have not yet been updated to run in parallel,
the multiprocessor check in init_atm_cases is still present.

* init_atmosphere/snow_albd_interp:
  Change scalefactor from c_float to RKIND
  Parallel interpolation of maximum snow albedo
  Enable the geotile manager to supersample tiles
  Initialize, finalize and push all modis snoalb onto the stack
  Update geotile manager to work with negative grid spacing
  Remove previous interpolation of snow albedo
  • Loading branch information
@mgduda
mgduda committed Dec 22, 2020
2 parents c8ca8fa + 9f4dfbf commit 057dd78
Show file tree
Hide file tree
Showing 2 changed files with 156 additions and 92 deletions.
34 changes: 25 additions & 9 deletions src/core_init_atmosphere/mpas_geotile_manager.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down
214 changes: 131 additions & 83 deletions src/core_init_atmosphere/mpas_init_atm_static.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -829,107 +832,152 @@ 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
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
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/'
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
do jTileStart = 1,02401,ny-6
jTileEnd = jTileStart + ny - 1 - 6
allocate(nhs(nCells))
allocate(snoalb_integer(nCells))
snoalb_integer(:) = 0
snoalb(:) = 0.0
nhs(:) = 0
fillval = 0.0
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)
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
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
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
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
do while (.not. mgr % is_stack_empty())
tile => mgr % pop_tile()
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
if (tile % is_processed) then
cycle
end if
iPoint = nearest_cell(lat_pt,lon_pt,iPoint,nCells,maxEdges, &
nEdgesOnCell,cellsOnCell, &
latCell,lonCell)
if (rarray(ii,jj,1) /= msgval) then
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
!
! 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
tile % is_processed = .true.
! 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
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
end do
end if
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
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(rarray)
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)', &
messageType=MPAS_LOG_CRIT)
endif
else if (trim(config_maxsnowalbedo_data) == 'NCEP') then
Expand Down

0 comments on commit 057dd78

Please sign in to comment.