Static interpolation psuedocode

static_interp.F

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! OG - Terrain static interpolation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! Create the KD Tree of nCells to be used for all static fields
allocate(kdcords(nCells))
do i = 1, nCells
   allocate(kdcords(i) % point(3))
   if (i > nCellsSolve) then
       kdcords(i) % is_halo = .true.
   else
       kdcords(i) % is_halo = .false.
enddo
tree => mpas_kd_construct(kdcords, 3)


! Call read the index file for this dataset
! create_pool(geog_pool)
call maps_parse_index(fname, geog_pool)

! From the goeg_pool, transfer the following variables:
call mpas_pool_get_config(geog_pool, 'dx', dx) ! The distance between two pixels
call mpas_pool_get_config(geog_pool, 'dy', dy) ! The distance between two pixels
call mpas_pool_get_config(geog_pool, 'tile_x', tile_nx) ! The number of pixels between two tiles
call mpas_pool_get_config(geog_pool, 'tile_y', tile_ny) ! The number of pixels between two tiles
call mpas_pool_get_config(geog_pool, 'tile_bdr', tile_bdr)
call mpas_pool_get_config(geog_pool, 'known_lat', start_lat)
call mpas_pool_get_config(geog_pool, 'known_lon', start_lon)
call mpas_pool_get_config(geog_pool, 'known_x', start_x)
call mpas_pool_get_config(geog_pool, 'known_y', start_y)
call mpas_pool_add_config(geog_pool, 'geog_dir_path', fname)

! Initalize the hash table by using dx, tile_nx, dy, tile_ny
nTileX = int((360.0_RKIND / dx)) / tile_nx
nTileY = int((180.0_RKIND / dy)) / tile_ny
allocate(hash(0:nTileX, 0:nTileY)) 

! Determine the amount of global pixels for this dataset
! (These are used in geotile_manager only.. so perhaps they should be
! there)
pixel_nx = 360.0_RKIND / dx
pixel_ny = 180.0_RKIND / dy
call mpas_pool_add_config(geog_pool, 'pixel_nx', pixel_nx)
call mpas_pool_add_config(geog_pool, 'pixel_ny', pixel_ny)

! Allocate the needed variables needed to interpolate the dataset
allocate(nhs(nCells))
allocate(ter(nCells))
nhs(:) = 0
ter(:) = 0
   
do cell = 1, nCells ! Loop through each cell to ensure each one has values

   if (nhs(cell) == 0) then
       ! If the tile hasn't been processed then get the tile and add it to the stck
       call mpas_get_tile(latCell(cell), lonCell(cell), tile, hash, geog_pool)
       stack => mpas_stack_push(stack, tile)
   endif

   do while ()
       tile => mpas_stack_pop_tile(stack)

       ! Tiles can be pushed onto the stack multiple times, don't process 
       ! it twice
       if ( tile % is_processed) then
           cycle
       else
           tile % is_processed = .true. ! Maybe move this after its actually proccesed.
       endif
  
      
       all_pixels_mapped_to_halo_cell = .true.


       do j in tile_ny + tile_bdr  
           do i in tile_ny + tile_bdr

               tval = tile % tile(i, j, 1)

               ! Calculate lat lon point of this datapoint
               lat_pt = real((j - (tile_bdr + 1) + tile % y(1) - 1), kind=RKIND) * dy + start_lat
               lon_pt = real((i - (tile_bdr + 1) + tile % x(1) - 1), kind=RKIND) * dx + start_lon = real(

               ! x,y,z = convert_lx(lat_pt, lon_pt)
               ! Search the kd tree for the cell that contains this point
               call mpas_kd_search(tree, (/x, y, z/), res)
               if (bdyMaskCell(res % cell) < nBdyLayers) then
                   ter(res % cell) = ter(res % cell) + tval
                   nhs(res % cell) = nhs(res % cell) + 1


                   if (res % cell <= nCellsSolve) then
                       all_pixels_mapped_to_halo_cell = .false.
                   endif
               else
                  ! Boundary cell code here
               endif
           enddo
       enddo
       if ( .not. all_pixels_mapped_to_halo_cell) then
           ! If all points didn't map to halo cells, then push the neighbor tiles 
           ! on the stack
           call mpas_get_tile_neighbors(tile, up, down, left, right, geog_pool)
           call push_tile_neighbors(stack, up, down, left, right, hash, geog_pool)
       endif
   enddo
enddo

do i = 1, nCells
    ter(i) = real(ter_integer(i)) / real(nhs(i))
enddo

deallocate(nhs)
deallocate(ter)

call mpas_stack_free(stack, free_payload=.true.)

! Use the hash table to deallocate each tile (if its associated)
do i = 0, nTileX
    do j = 0, nTileY
        if (assocaited(hash(i,j) % ptr)) then
            deallocate(hash(i,j) % ptr % tile)
            deallocate(hash(i,j) % ptr)
        endif
    enddo
enddo
deallocate(hash)

call mpas_pool_destroy_pool(geog_pool)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! New - Terrain static interpolation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! Create the KD Tree of nCells to be used for all static fields
allocate(kdcords(nCells))
do i = 1, nCells
   allocate(kdcords(i) % point(3))
   if (i > nCellsSolve) then
       kdcords(i) % is_halo = .true.
   else
       kdcords(i) % is_halo = .false.
enddo
tree => mpas_kd_construct(kdcords, 3)

! Get the directory for the geodata set from the
! namelist
index_file_name = namelist(directory_name)

! Initalize geotile manager for this geodata set
! geotile_manager_init
!    * call mpas_parse_geoindex and create pool
!    * Calculate values needed by the geotile manager (hash table size)
!       - Hash table size
!       - Values needed to convert pixels cords to lat lon, and vice versa
!    * Allocate the hash table
ierr = mpas_geotile_manager_init(index_file_name)

! Somehow keep track of the state of the geotile manager, either return
! a context-like object, or another method?

! Somehow we will need to get the following variables from the geotile
! manager:
! - tile_nx, tile_ny
! - tile_bdr
! - dx, dy
! - start_lat, start_lon
! 
! I think we could have the geotile manager return these fields in a data
! structure, but I think it makes most sense to return the pool thats created
! by mpas_parse_geoindex. That way, if we, or perhaps another use needs another
! variable thats defined in the index file, they can retrive it.
!

allocate(nhs(nCells))
allocate(ter(nCells))

nhs(:) = 0
ter(:) = 0

do cell = 1, nCells
    if (nhs(cell) == 0) then
       call mpas_get_tile(latCell(cell), lonCell(cell), tile)
       stack => mpas_stack_push(stack, tile)
    endif

    do while(.not. mpas_stack_is_empty(stack))
       if (tile % is_processed) then
           cycle
       else
           tile % is_processed = .true.
       endif

       all_pixels_mapped_to_halo_cell = .true.

       do j = 1 + tile_bdr, tile_ny + tile_bdr
           do i = 1 + tile_bdr, tile_nx + tile_bdr
               tval = tile % tile(i, j, 1)

               lat_pt = real((j - (tile_bdr + 1) + tile % y(1) - 1), kind=RKIND) * dy + start_lat
               lon_pt = real((i - (tile_bdr + 1) + tile % x(1) - 1), kind=RKIND) * dx + start_lon
               
               ! convert to radians

               call convert_lx(x, y, z, sphere_radius, lat_pt, lon_pt)
               call mpas_kd_search(tree, (/x, y, z/), res)
               
               if (bdyMaskCell(res % cell) < nBdyLayers) then
                   ter(res % cell) = ter(res % cell) + tval
                   nhs(res % cell) = nhs (res % cell) + 1

                   if (res % cell <= nCellsSolve) then
                       all_pixels_mapped_to_halo_cell = .false.
                   endif
               else
                  ! Code for boundary cells here
               endif
           enddo
       enddo

       ! Push the tile neighbors on the cell, if we need to
       if (.not. all_pixels_mapped_to_halo_cell) then
           call mpas_get_tile_neighbors(tile, up, down, left, right)
           call push_tile_neighbors(stack, up, down, left, right)
       endif
    enddo
enddo

do i = 1, nCells
   ter(i) = real(ter_integer(i)) / real(nhs(i))
enddo

deallocate(nhs)
deallocate(ter)

call mpas_stack_free(stack)
call mpas_geotile_manager_free() ! Argument ???