Merge branch 'init_atmosphere/snow_albd_interp' into develop (PR #750)

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

Author: mgduda

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
@@ -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)

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
 
@@ -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
 
@@ -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