/WPS/geogrid/src/process_tile_module.F

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Module: process_tile
!
! Description:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module process_tile_module

   use module_debug


   contains


   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: process_tile
   !
   ! Purpose: To process a tile, whose lower-left corner is at 
   !       (tile_i_min, tile_j_min) and whose upper-right corner is at
   !       (tile_i_max, tile_j_max), of the model grid given by which_domain 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   subroutine process_tile(which_domain, grid_type, dynopt,        &
                           dummy_start_dom_i, dummy_end_dom_i,     &
                           dummy_start_dom_j, dummy_end_dom_j,     &
                           dummy_start_patch_i, dummy_end_patch_i, &
                           dummy_start_patch_j, dummy_end_patch_j, &
                           extra_col, extra_row)
   
      use bitarray_module
      use hash_module
      use llxy_module
      use misc_definitions_module
      use output_module
      use smooth_module
      use source_data_module
      use gridinfo_module, only : subgrid_ratio_x, subgrid_ratio_y
    
      implicit none

      ! Arguments
      integer, intent(in) :: which_domain, dynopt, &
                             dummy_start_dom_i, dummy_end_dom_i, dummy_start_dom_j, dummy_end_dom_j, &
                             dummy_start_patch_i, dummy_end_patch_i, dummy_start_patch_j, dummy_end_patch_j
      logical, intent(in) :: extra_col, extra_row
      character (len=1), intent(in) :: grid_type
    
      ! Local variables
      integer :: i, j, k, kk, istatus, ifieldstatus, idomcatstatus, field_count
      integer :: min_category, max_category, min_level, max_level, &
                 smth_opt, smth_passes, num_landmask_categories
      integer :: start_dom_i, end_dom_i, start_dom_j, end_dom_j, end_dom_stag_i, end_dom_stag_j
      integer :: start_patch_i, end_patch_i, start_patch_j, end_patch_j, end_patch_stag_i, end_patch_stag_j
      integer :: start_mem_i, end_mem_i, start_mem_j, end_mem_j, end_mem_stag_i, end_mem_stag_j
      integer :: start_mem_sub_i, end_mem_sub_i, start_mem_sub_j, end_mem_sub_j
      integer :: sm1, em1, sm2, em2
      integer :: istagger
      integer, dimension(MAX_LANDMASK_CATEGORIES) :: landmask_value
      real :: sum, dominant, msg_fill_val, topo_flag_val, mass_flag, land_total, water_total
      real, dimension(16) :: corner_lats, corner_lons
      real, pointer, dimension(:,:) :: xlat_array,   xlon_array,   &
                                       xlat_array_u, xlon_array_u, &
                                       xlat_array_v, xlon_array_v, &
                                       clat_array,   clon_array,   &
                                       xlat_array_subgrid, xlon_array_subgrid, &
                                       f_array, e_array, &
                                       mapfac_array_m_x, mapfac_array_u_x, mapfac_array_v_x, &
                                       mapfac_array_m_y, mapfac_array_u_y, mapfac_array_v_y, &
                                       mapfac_array_x_subgrid, mapfac_array_y_subgrid,       &
                                       sina_array, cosa_array
      real, pointer, dimension(:,:) :: xlat_ptr, xlon_ptr, mapfac_ptr_x, mapfac_ptr_y, landmask, dominant_field
      real, pointer, dimension(:,:,:) :: field, slp_field
      logical :: is_water_mask, only_save_dominant, halt_on_missing
      character (len=19) :: datestr
      character (len=128) :: fieldname, gradname, domname, landmask_name
      character (len=256) :: temp_string
      type (bitarray) :: processed_domain 
      type (hashtable) :: processed_fieldnames
      character (len=128), dimension(2) :: dimnames
      integer :: sr_x, sr_y
      logical :: subgrid_var
  
      sr_x=subgrid_ratio_x(which_domain)
      sr_y=subgrid_ratio_y(which_domain)

      datestr = '0000-00-00_00:00:00'
      field_count = 0
      mass_flag=1.0
    
      ! The following pertains primarily to the C grid
      ! Determine whether only (n-1)th rows/columns should be computed for variables
      !   on staggered grid. In a distributed memory situation, not every tile should
      !   have only (n-1)th rows/columns computed, or we end up with (n-k) 
      !   rows/columns when there are k patches in the y/x direction
      if (extra_col) then
         start_patch_i    = dummy_start_patch_i    ! The seemingly pointless renaming of start
         end_patch_i      = dummy_end_patch_i - 1  !   naming convention with modified end_patch variables, 
         end_patch_stag_i = dummy_end_patch_i      !   variables is so that we can maintain consistent
                                                   !   which are marked as intent(in)
         start_mem_i    = start_patch_i    - HALO_WIDTH
         end_mem_i      = end_patch_i      + HALO_WIDTH
         end_mem_stag_i = end_patch_stag_i + HALO_WIDTH
      else                                     
         start_patch_i    = dummy_start_patch_i
         end_patch_i      = dummy_end_patch_i
         end_patch_stag_i = dummy_end_patch_i

         start_mem_i    = start_patch_i    - HALO_WIDTH
         end_mem_i      = end_patch_i      + HALO_WIDTH
         end_mem_stag_i = end_patch_stag_i + HALO_WIDTH
      end if
    
      if (extra_row) then
         start_patch_j    = dummy_start_patch_j
         end_patch_j      = dummy_end_patch_j - 1
         end_patch_stag_j = dummy_end_patch_j

         start_mem_j    = start_patch_j    - HALO_WIDTH
         end_mem_j      = end_patch_j      + HALO_WIDTH
         end_mem_stag_j = end_patch_stag_j + HALO_WIDTH
      else
         start_patch_j    = dummy_start_patch_j
         end_patch_j      = dummy_end_patch_j
         end_patch_stag_j = dummy_end_patch_j

         start_mem_j    = start_patch_j    - HALO_WIDTH
         end_mem_j      = end_patch_j      + HALO_WIDTH
         end_mem_stag_j = end_patch_stag_j + HALO_WIDTH
      end if

      start_dom_i = dummy_start_dom_i
      if (grid_type == 'C') then
         end_dom_i      = dummy_end_dom_i - 1
         end_dom_stag_i = dummy_end_dom_i
      else if (grid_type == 'E') then
         end_dom_i      = dummy_end_dom_i
         end_dom_stag_i = dummy_end_dom_i
      end if

      start_dom_j = dummy_start_dom_j
      if (grid_type == 'C') then
         end_dom_j      = dummy_end_dom_j - 1
         end_dom_stag_j = dummy_end_dom_j
      else if (grid_type == 'E') then
         end_dom_j      = dummy_end_dom_j
         end_dom_stag_j = dummy_end_dom_j
      end if

      start_mem_sub_i = (start_mem_i - 1 ) * sr_x + 1
      if(.not.extra_col)then
        end_mem_sub_i = (end_mem_i) * sr_x
      else
        end_mem_sub_i = (end_mem_i + 1) * sr_x
      end if

      start_mem_sub_j = (start_mem_j -1 ) * sr_y + 1
      if(.not.extra_row)then
        end_mem_sub_j = (end_mem_j) * sr_y
      else
        end_mem_sub_j = (end_mem_j + 1) * sr_y
      end if
    
      ! Allocate arrays to hold all lat/lon fields; these will persist for the duration of
      !   the process_tile routine
      ! For C grid, we have M, U, and V points
      ! For E grid, we have only M and V points
      allocate(xlat_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
      allocate(xlon_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
      allocate(xlat_array_v(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
      allocate(xlon_array_v(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
      if (grid_type == 'C') then
         allocate(xlat_array_u(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
         allocate(xlon_array_u(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
         allocate(clat_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
         allocate(clon_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
         allocate(xlat_array_subgrid(start_mem_sub_i:end_mem_sub_i,start_mem_sub_j:end_mem_sub_j))
         allocate(xlon_array_subgrid(start_mem_sub_i:end_mem_sub_i,start_mem_sub_j:end_mem_sub_j))
         allocate(mapfac_array_x_subgrid(start_mem_sub_i:end_mem_sub_i,start_mem_sub_j:end_mem_sub_j))
         allocate(mapfac_array_y_subgrid(start_mem_sub_i:end_mem_sub_i,start_mem_sub_j:end_mem_sub_j))
      end if
    
      ! Initialize hash table to track which fields have been processed
      call hash_init(processed_fieldnames)
    
      !
      ! Calculate lat/lon for every point in the tile (XLAT and XLON)
      ! The xlat_array and xlon_array arrays will be used in processing other fields
      !
      call mprintf(.true.,STDOUT,'  Processing XLAT and XLONG')
    
      if (grid_type == 'C') then
         call get_lat_lon_fields(xlat_array, xlon_array, start_mem_i, &
                               start_mem_j, end_mem_i, end_mem_j, M)
         call get_lat_lon_fields(xlat_array_v, xlon_array_v, start_mem_i, &
                               start_mem_j, end_mem_i, end_mem_stag_j, V)
         call get_lat_lon_fields(xlat_array_u, xlon_array_u, start_mem_i, &
                               start_mem_j, end_mem_stag_i, end_mem_j, U)
         call get_lat_lon_fields(clat_array, clon_array, start_mem_i, &
                               start_mem_j, end_mem_i, end_mem_j, M, comp_ll=.true.)
         call get_lat_lon_fields(xlat_array_subgrid, xlon_array_subgrid, &
                                 start_mem_sub_i, start_mem_sub_j, end_mem_sub_i, end_mem_sub_j, M, sub_x=sr_x, sub_y=sr_y)
         call get_map_factor(xlat_array_subgrid, xlon_array_subgrid, mapfac_array_x_subgrid, &
                             mapfac_array_y_subgrid, start_mem_sub_i, start_mem_sub_j, end_mem_sub_i, end_mem_sub_j)

         corner_lats(1) = xlat_array(start_patch_i,start_patch_j)
         corner_lats(2) = xlat_array(start_patch_i,end_patch_j)
         corner_lats(3) = xlat_array(end_patch_i,end_patch_j)
         corner_lats(4) = xlat_array(end_patch_i,start_patch_j)
     
         corner_lats(5) = xlat_array_u(start_patch_i,start_patch_j)
         corner_lats(6) = xlat_array_u(start_patch_i,end_patch_j)
         corner_lats(7) = xlat_array_u(end_patch_stag_i,end_patch_j)
         corner_lats(8) = xlat_array_u(end_patch_stag_i,start_patch_j)
     
         corner_lats(9)  = xlat_array_v(start_patch_i,start_patch_j)
         corner_lats(10) = xlat_array_v(start_patch_i,end_patch_stag_j)
         corner_lats(11) = xlat_array_v(end_patch_i,end_patch_stag_j)
         corner_lats(12) = xlat_array_v(end_patch_i,start_patch_j)
     
         call xytoll(real(start_patch_i)-0.5, real(start_patch_j)-0.5, corner_lats(13), corner_lons(13), M)
         call xytoll(real(start_patch_i)-0.5, real(end_patch_j)+0.5, corner_lats(14), corner_lons(14), M)
         call xytoll(real(end_patch_i)+0.5, real(end_patch_j)+0.5, corner_lats(15), corner_lons(15), M)
         call xytoll(real(end_patch_i)+0.5, real(start_patch_j)-0.5, corner_lats(16), corner_lons(16), M)

         corner_lons(1) = xlon_array(start_patch_i,start_patch_j)
         corner_lons(2) = xlon_array(start_patch_i,end_patch_j)
         corner_lons(3) = xlon_array(end_patch_i,end_patch_j)
         corner_lons(4) = xlon_array(end_patch_i,start_patch_j)
     
         corner_lons(5) = xlon_array_u(start_patch_i,start_patch_j)
         corner_lons(6) = xlon_array_u(start_patch_i,end_patch_j)
         corner_lons(7) = xlon_array_u(end_patch_stag_i,end_patch_j)
         corner_lons(8) = xlon_array_u(end_patch_stag_i,start_patch_j)
     
         corner_lons(9)  = xlon_array_v(start_patch_i,start_patch_j)
         corner_lons(10) = xlon_array_v(start_patch_i,end_patch_stag_j)
         corner_lons(11) = xlon_array_v(end_patch_i,end_patch_stag_j)
         corner_lons(12) = xlon_array_v(end_patch_i,start_patch_j)
     
      else if (grid_type == 'E') then
         call get_lat_lon_fields(xlat_array, xlon_array, start_mem_i, &
                               start_mem_j, end_mem_i, end_mem_j, HH)
         call get_lat_lon_fields(xlat_array_v, xlon_array_v, start_mem_i, &
                               start_mem_j, end_mem_i, end_mem_stag_j, VV)
   
         corner_lats(1) = xlat_array(start_patch_i,start_patch_j)
         corner_lats(2) = xlat_array(start_patch_i,end_patch_j)
         corner_lats(3) = xlat_array(end_patch_i,end_patch_j)
         corner_lats(4) = xlat_array(end_patch_i,start_patch_j)
     
         corner_lats(5) = xlat_array_v(start_patch_i,start_patch_j)
         corner_lats(6) = xlat_array_v(start_patch_i,end_patch_stag_j)
         corner_lats(7) = xlat_array_v(end_patch_i,end_patch_stag_j)
         corner_lats(8) = xlat_array_v(end_patch_i,start_patch_j)
     
         corner_lats(9)  = 0.0
         corner_lats(10) = 0.0
         corner_lats(11) = 0.0
         corner_lats(12) = 0.0
     
         corner_lats(13) = 0.0
         corner_lats(14) = 0.0
         corner_lats(15) = 0.0
         corner_lats(16) = 0.0
     
         corner_lons(1) = xlon_array(start_patch_i,start_patch_j)
         corner_lons(2) = xlon_array(start_patch_i,end_patch_j)
         corner_lons(3) = xlon_array(end_patch_i,end_patch_j)
         corner_lons(4) = xlon_array(end_patch_i,start_patch_j)
     
         corner_lons(5) = xlon_array_v(start_patch_i,start_patch_j)
         corner_lons(6) = xlon_array_v(start_patch_i,end_patch_stag_j)
         corner_lons(7) = xlon_array_v(end_patch_i,end_patch_stag_j)
         corner_lons(8) = xlon_array_v(end_patch_i,start_patch_j)
     
         corner_lons(9)  = 0.0
         corner_lons(10) = 0.0
         corner_lons(11) = 0.0
         corner_lons(12) = 0.0
     
         corner_lons(13) = 0.0
         corner_lons(14) = 0.0
         corner_lons(15) = 0.0
         corner_lons(16) = 0.0
    
      end if

      ! Initialize the output module now that we have the corner point lats/lons
      call output_init(which_domain, 'OUTPUT FROM GEOGRID V3.3', '0000-00-00_00:00:00', grid_type, dynopt, &
                       corner_lats, corner_lons, &
                       start_dom_i,   end_dom_i,   start_dom_j,   end_dom_j, &
                       start_patch_i, end_patch_i, start_patch_j, end_patch_j, &
                       start_mem_i,   end_mem_i,   start_mem_j,   end_mem_j, &
                       extra_col, extra_row)
    
      call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
                       'XLAT_M', datestr, real_array = xlat_array)
      call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
                       'XLONG_M', datestr, real_array = xlon_array)
      call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, &
                       'XLAT_V', datestr, real_array = xlat_array_v)
      call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, &
                       'XLONG_V', datestr, real_array = xlon_array_v)
      if (grid_type == 'C') then
         call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, &
                        'XLAT_U', datestr, real_array = xlat_array_u)
         call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, &
                        'XLONG_U', datestr, real_array = xlon_array_u)
         call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
                        'CLAT', datestr, real_array = clat_array)
         call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
                        'CLONG', datestr, real_array = clon_array)
         call write_field(start_mem_sub_i, end_mem_sub_i, start_mem_sub_j, end_mem_sub_j, 1, 1, &
                        'FXLAT', datestr, real_array = xlat_array_subgrid)
         call write_field(start_mem_sub_i, end_mem_sub_i, start_mem_sub_j, end_mem_sub_j, 1, 1, &
                        'FXLONG', datestr, real_array = xlon_array_subgrid)

         if (associated(clat_array)) deallocate(clat_array)
         if (associated(clon_array)) deallocate(clon_array)

      end if
   

      !
      ! Calculate map factor for current domain
      !
      if (grid_type == 'C') then
         call mprintf(.true.,STDOUT,'  Processing MAPFAC')
    
         allocate(mapfac_array_m_x(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
         allocate(mapfac_array_m_y(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
         call get_map_factor(xlat_array, xlon_array, mapfac_array_m_x, mapfac_array_m_y, start_mem_i, &
                             start_mem_j, end_mem_i, end_mem_j)
! Global WRF uses map scale factors in X and Y directions, but "regular" WRF uses a single MSF
!    on each staggering. In the case of regular WRF, we can assume that MAPFAC_MX = MAPFAC_MY = MAPFAC_M,
!    and so we can simply write MAPFAC_MX as the MAPFAC_M field. Ultimately, when global WRF is
!    merged into the WRF trunk, we will need only two map scale factor fields for each staggering,
!    in the x and y directions, and these will be the same in the case of non-Cassini projections
         call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_M',  &
                          datestr, real_array = mapfac_array_m_x)
         call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_MX', &
                          datestr, real_array = mapfac_array_m_x)
         call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_MY', &
                         datestr, real_array = mapfac_array_m_y)
    
         allocate(mapfac_array_v_x(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
         allocate(mapfac_array_v_y(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
         call get_map_factor(xlat_array_v, xlon_array_v, mapfac_array_v_x, mapfac_array_v_y, start_mem_i, &
                             start_mem_j, end_mem_i, end_mem_stag_j)
         call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'MAPFAC_V',  &
                          datestr, real_array = mapfac_array_v_x)
         call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'MAPFAC_VX', &
                          datestr, real_array = mapfac_array_v_x)
         call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'MAPFAC_VY', &
                          datestr, real_array = mapfac_array_v_y)
    
         allocate(mapfac_array_u_x(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
         allocate(mapfac_array_u_y(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
         call get_map_factor(xlat_array_u, xlon_array_u, mapfac_array_u_x, mapfac_array_u_y, start_mem_i, &
                             start_mem_j, end_mem_stag_i, end_mem_j)
         call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_U',  &
                          datestr, real_array = mapfac_array_u_x)
         call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_UX', &
                          datestr, real_array = mapfac_array_u_x)
         call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_UY', &
                          datestr, real_array = mapfac_array_u_y)

      end if
    
    
      !
      ! Coriolis parameters (E and F)
      !
      call mprintf(.true.,STDOUT,'  Processing F and E')
    
      allocate(f_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
      allocate(e_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
    
      call get_coriolis_parameters(xlat_array, f_array, e_array, &
                                   start_mem_i, start_mem_j, end_mem_i, end_mem_j)
    
      call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'E', &
                       datestr, real_array = e_array)
      call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'F', &
                       datestr, real_array = f_array)
    
      if (associated(f_array)) deallocate(f_array)
      if (associated(e_array)) deallocate(e_array)
    
    
      !
      ! Rotation angle (SINALPHA and COSALPHA)
      !
      if (grid_type == 'C') then
         call mprintf(.true.,STDOUT,'  Processing ROTANG')
    
         allocate(sina_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
         allocate(cosa_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
    
         call get_rotang(xlat_array, xlon_array, cosa_array, sina_array, &
                         start_mem_i, start_mem_j, end_mem_i, end_mem_j)
    
         call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'SINALPHA', &
                          datestr, real_array = sina_array)
         call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'COSALPHA', &
                          datestr, real_array = cosa_array)
    
         if (associated(sina_array)) deallocate(sina_array)
         if (associated(cosa_array)) deallocate(cosa_array)
      end if
    
      ! Every field up until now should probably just be processed regardless of what the user 
      !   has specified for fields to be processed.
      ! Hereafter, we process user-specified fields
    
      !
      ! First process the field that we will derive a landmask from
      !
      call get_landmask_field(geog_data_res(which_domain), landmask_name, is_water_mask, landmask_value, istatus)

      do kk=1,MAX_LANDMASK_CATEGORIES
         if (landmask_value(kk) == INVALID) then
            num_landmask_categories = kk-1
            exit
         end if
      end do
      if (kk > MAX_LANDMASK_CATEGORIES) num_landmask_categories = MAX_LANDMASK_CATEGORIES

      if (istatus /= 0) then
         call mprintf(.true.,WARN,'No field specified for landmask calculation. Will set landmask=1 at every grid point.')
     
         allocate(landmask(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
         landmask = 1.
         call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'LANDMASK', &
                          datestr, landmask)
    
      else
    
         allocate(landmask(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
         landmask = 1.
     
         call mprintf(.true.,STDOUT,'  Processing %s', s1=trim(landmask_name))
     
         call get_missing_fill_value(landmask_name, msg_fill_val, istatus)
         if (istatus /= 0) msg_fill_val = NAN
     
         call get_halt_on_missing(landmask_name, halt_on_missing, istatus)
         if (istatus /= 0) halt_on_missing = .false.
     
         ! Do we calculate a dominant category for this field?
         call get_domcategory_name(landmask_name, domname, only_save_dominant, idomcatstatus)
     
         temp_string = ' '
         temp_string(1:128) = landmask_name
         call hash_insert(processed_fieldnames, temp_string)
     
         call get_max_categories(landmask_name, min_category, max_category, istatus)
         allocate(field(start_mem_i:end_mem_i, start_mem_j:end_mem_j, min_category:max_category))

         if (.not. only_save_dominant) then
            field_count = field_count + 1
            call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
                         i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=landmask_name)
         else
            field_count = field_count + 1
            call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
                         i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
         end if

         if (grid_type == 'C') then
            call calc_field(landmask_name, field, xlat_array, xlon_array, M, &
                            start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
                            min_category, max_category, processed_domain, 1, landmask=landmask, sr_x=1, sr_y=1)
         else if (grid_type == 'E') then
            call calc_field(landmask_name, field, xlat_array, xlon_array, HH, &
                            start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
                            min_category, max_category, processed_domain, 1, landmask=landmask, sr_x=1, sr_y=1)
         end if
     
         ! Find fractions
         do i=start_mem_i, end_mem_i
            do j=start_mem_j, end_mem_j
               sum = 0.
               do k=min_category,max_category
                  sum = sum + field(i,j,k)
               end do
               if (sum > 0.0) then
                  do k=min_category,max_category
                     field(i,j,k) = field(i,j,k) / sum
                  end do
               else
                  do k=min_category,max_category
                     field(i,j,k) = msg_fill_val
                  end do
               end if
            end do
         end do
     
         ! If user wants to halt when a missing value is found in output field, check now
         if (halt_on_missing) then
            do i=start_mem_i, end_mem_i
               do j=start_mem_j, end_mem_j
                  ! Only need to examine k=1
                  if (field(i,j,1) == msg_fill_val) then
                     call mprintf(.true.,ERROR,' Missing value encountered in output field. Quitting.')
                  end if
               end do
            end do
         end if
     
         if (is_water_mask) then
            call mprintf(.true.,STDOUT,'  Calculating landmask from %s ( WATER =', &
                         newline=.false.,s1=trim(landmask_name))
         else
            call mprintf(.true.,STDOUT,'  Calculating landmask from %s ( LAND =', &
                         newline=.false.,s1=trim(landmask_name))
         end if
         do k = 1, num_landmask_categories
            call mprintf(.true.,STDOUT,' %i',newline=.false.,i1=landmask_value(k))
            if (k == num_landmask_categories) call mprintf(.true.,STDOUT,')')
         end do
     
         ! Calculate landmask
         if (is_water_mask) then
            do i=start_mem_i, end_mem_i
               do j=start_mem_j, end_mem_j
                  water_total = -1.
                  do k=1,num_landmask_categories
                     if (landmask_value(k) >= min_category .and. landmask_value(k) <= max_category) then 
                        if (field(i,j,landmask_value(k)) /= msg_fill_val) then
                           if (water_total < 0.) water_total = 0.
                           water_total = water_total + field(i,j,landmask_value(k)) 
                        else
                           water_total = -1.
                           exit
                        end if
                     end if
                  end do
                  if (water_total >= 0.0) then
                     if (water_total < 0.50) then
                        landmask(i,j) = 1.
                     else
                        landmask(i,j) = 0.
                     end if
                  else
                     landmask(i,j) = -1.
                  end if
               end do
            end do
         else
            do i=start_mem_i, end_mem_i
               do j=start_mem_j, end_mem_j
                  land_total = -1.
                  do k=1,num_landmask_categories
                     if (landmask_value(k) >= min_category .and. landmask_value(k) <= max_category) then 
                        if (field(i,j,landmask_value(k)) /= msg_fill_val) then
                           if (land_total < 0.) land_total = 0.
                           land_total = land_total + field(i,j,landmask_value(k)) 
                        else
                           land_total = -1.
                           exit
                        end if
                     end if
                  end do
                  if (land_total >= 0.0) then
                     if (land_total > 0.50) then
                        landmask(i,j) = 1.
                     else
                        landmask(i,j) = 0.
                     end if
                  else
                     landmask(i,j) = -1.
                  end if
               end do
            end do
         end if
    
         call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'LANDMASK', &
                          datestr, landmask)
     
         ! If we should only save the dominant category, then no need to write out fractional field
         if (.not.only_save_dominant .or. (idomcatstatus /= 0)) then
     
            ! Finally, we may be asked to smooth the fractional field
            call get_smooth_option(landmask_name, smth_opt, smth_passes, istatus)
            if (istatus == 0) then

               if (grid_type == 'C') then
                  if (smth_opt == ONETWOONE) then
                     call one_two_one(field,                      &
                                      start_patch_i, end_patch_i, &
                                      start_patch_j, end_patch_j, &
                                      start_mem_i, end_mem_i,     &
                                      start_mem_j, end_mem_j,     &
                                      min_category, max_category, &
                                      smth_passes, msg_fill_val)
                  else if (smth_opt == SMTHDESMTH) then
                     call smth_desmth(field,                      &
                                      start_patch_i, end_patch_i, &
                                      start_patch_j, end_patch_j, &
                                      start_mem_i, end_mem_i,     &
                                      start_mem_j, end_mem_j,     &
                                      min_category, max_category, &
                                      smth_passes, msg_fill_val)
                  else if (smth_opt == SMTHDESMTH_SPECIAL) then
                     call smth_desmth_special(field,              &
                                      start_patch_i, end_patch_i, &
                                      start_patch_j, end_patch_j, &
                                      start_mem_i, end_mem_i,     &
                                      start_mem_j, end_mem_j,     &
                                      min_category, max_category, &
                                      smth_passes, msg_fill_val)
                  end if
               else if (grid_type == 'E') then
                  if (smth_opt == ONETWOONE) then
                     call one_two_one_egrid(field,                &
                                      start_patch_i, end_patch_i, &
                                      start_patch_j, end_patch_j, &
                                      start_mem_i, end_mem_i,     &
                                      start_mem_j, end_mem_j,     &
                                      min_category, max_category, &
                                      smth_passes, msg_fill_val, 1.0)
                  else if (smth_opt == SMTHDESMTH) then
                     call smth_desmth_egrid(field,                &
                                      start_patch_i, end_patch_i, &
                                      start_patch_j, end_patch_j, &
                                      start_mem_i, end_mem_i,     &
                                      start_mem_j, end_mem_j,     &
                                      min_category, max_category, &
                                      smth_passes, msg_fill_val, 1.0)
                  else if (smth_opt == SMTHDESMTH_SPECIAL) then
                     call mprintf(.true.,WARN,'smth-desmth_special is not currently implemented for NMM. '// &
                                              'No smoothing will be done.')
                  end if
               end if

            end if
      
            call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
                             min_category, max_category, trim(landmask_name), &
                             datestr, real_array=field)
         end if
     
         if (idomcatstatus == 0) then
            allocate(dominant_field(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
     
            if (.not. only_save_dominant) then
               field_count = field_count + 1
               call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
                            i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
            end if

            do i=start_mem_i, end_mem_i
               do j=start_mem_j, end_mem_j
                  if ((landmask(i,j) == 1. .and. is_water_mask) .or. &
                      (landmask(i,j) == 0. .and. .not.is_water_mask)) then
                     dominant = 0.
                     dominant_field(i,j) = real(min_category-1)
                     do k=min_category,max_category
                        do kk=1,num_landmask_categories
                           if (k == landmask_value(kk)) exit 
                        end do
                        if (field(i,j,k) > dominant .and. kk > num_landmask_categories) then
                           dominant_field(i,j) = real(k)
                           dominant = field(i,j,k)
                        end if
                     end do
                  else
                     dominant = 0.
                     dominant_field(i,j) = real(min_category-1)
                     do k=min_category,max_category
                        do kk=1,num_landmask_categories
                           if (field(i,j,k) > dominant .and. k == landmask_value(kk)) then
                              dominant_field(i,j) = real(k)
                              dominant = field(i,j,k)
                           end if
                        end do
                     end do
                  end if
               end do
            end do
     
            call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, trim(domname), &
                             datestr, dominant_field)
          
            deallocate(dominant_field)
         end if
     
         deallocate(field)
      end if
   
      !
      ! Now process all other fields specified by the user
      !
      call reset_next_field()
      ifieldstatus = 0
      do while (ifieldstatus == 0) 
         call get_next_fieldname(fieldname, ifieldstatus)
     
         ! There is another field in the GEOGRID.TBL file 
         if (ifieldstatus == 0) then
            temp_string(1:128) = fieldname
      
            ! If this field is still to be processed
            if (.not. hash_search(processed_fieldnames, temp_string)) then
     
               call hash_insert(processed_fieldnames, temp_string)
               call mprintf(.true.,STDOUT,'  Processing %s', s1=trim(fieldname))
       
               call get_output_stagger(fieldname, istagger, istatus)
               dimnames(:) = 'null'

               subgrid_var=is_subgrid_var(fieldname)
               if( subgrid_var ) call get_subgrid_dim_name(dimnames)
       
               if (istagger == M .and. .not. subgrid_var ) then
                  sm1 = start_mem_i
                  em1 = end_mem_i
                  sm2 = start_mem_j
                  em2 = end_mem_j
                  xlat_ptr => xlat_array 
                  xlon_ptr => xlon_array 
                  mapfac_ptr_x => mapfac_array_m_x
                  mapfac_ptr_y => mapfac_array_m_y
               else if ( subgrid_var ) then
                  sm1 = start_mem_sub_i
                  em1 = end_mem_sub_i
                  sm2 = start_mem_sub_j
                  em2 = end_mem_sub_j
                  xlat_ptr => xlat_array_subgrid
                  xlon_ptr => xlon_array_subgrid
                  mapfac_ptr_x => mapfac_array_x_subgrid
                  mapfac_ptr_y => mapfac_array_y_subgrid
               else if (istagger == U) then ! In the case that extra_cols = .false.
                  sm1 = start_mem_i          ! we should have that end_mem_stag is
                  em1 = end_mem_stag_i       ! the same as end_mem, so we do not need
                  sm2 = start_mem_j          ! to check extra_cols or extra rows here
                  em2 = end_mem_j
                  xlat_ptr => xlat_array_u 
                  xlon_ptr => xlon_array_u 
                  mapfac_ptr_x => mapfac_array_u_x
                  mapfac_ptr_y => mapfac_array_u_y
               else if (istagger == V) then
                  sm1 = start_mem_i
                  em1 = end_mem_i
                  sm2 = start_mem_j
                  em2 = end_mem_stag_j
                  xlat_ptr => xlat_array_v 
                  xlon_ptr => xlon_array_v 
                  mapfac_ptr_x => mapfac_array_v_x
                  mapfac_ptr_y => mapfac_array_v_y
               else if (istagger == HH) then   ! E grid
                  sm1 = start_mem_i
                  em1 = end_mem_i
                  sm2 = start_mem_j
                  em2 = end_mem_j
                  xlat_ptr => xlat_array 
                  xlon_ptr => xlon_array 
                  mapfac_ptr_x => mapfac_array_m_x
                  mapfac_ptr_y => mapfac_array_m_y
               else if (istagger == VV) then   ! E grid 
                  sm1 = start_mem_i
                  em1 = end_mem_i
                  sm2 = start_mem_j
                  em2 = end_mem_stag_j
                  xlat_ptr => xlat_array_v 
                  xlon_ptr => xlon_array_v 
                  mapfac_ptr_x => mapfac_array_v_x
                  mapfac_ptr_y => mapfac_array_v_y
               end if

      
               call get_missing_fill_value(fieldname, msg_fill_val, istatus)
               if (istatus /= 0) msg_fill_val = NAN 
       
               call get_halt_on_missing(fieldname, halt_on_missing, istatus)
               if (istatus /= 0) halt_on_missing = .false.
       
               ! Destination field type is CONTINUOUS
               if (iget_fieldtype(fieldname,istatus) == CONTINUOUS) then
                  call get_max_levels(fieldname, min_level, max_level, istatus)
                  allocate(field(sm1:em1, sm2:em2, min_level:max_level))

                  field_count = field_count + 1
                  call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
                               i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=fieldname)

                  if ( subgrid_var ) then
                     call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
                                sm1, em1, sm2, em2, min_level, max_level, &
                                processed_domain, 1, sr_x=sr_x, sr_y=sr_y)
                  else
                     call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
                                sm1, em1, sm2, em2, min_level, max_level, &
                                processed_domain, 1, landmask=landmask)
                  end if
        
                  ! If user wants to halt when a missing value is found in output field, check now
                  if (halt_on_missing) then
                     do i=sm1, em1
                        do j=sm2, em2
                           ! Only need to examine k=1
                           if (field(i,j,1) == msg_fill_val) then
                              call mprintf(.true.,ERROR,' Missing value encountered in output field. Quitting.')
                           end if
                        end do
                     end do
                  end if
        
                  ! We may be asked to smooth the fractional field
                  call get_smooth_option(fieldname, smth_opt, smth_passes, istatus)
                  if (istatus == 0) then

                     if (grid_type == 'C') then
                        if (smth_opt == ONETWOONE) then
                           call one_two_one(field,                      &
                                            start_patch_i, end_patch_i, &
                                            start_patch_j, end_patch_j, &
                                            sm1,         em1,           &
                                            sm2,         em2,           &
                                            min_level, max_level,       &
                                            smth_passes, msg_fill_val)
                        else if (smth_opt == SMTHDESMTH) then
                           call smth_desmth(field,                      &
                                            start_patch_i, end_patch_i, &
                                            start_patch_j, end_patch_j, &
                                            sm1,         em1,           &
                                            sm2,         em2,           &
                                            min_level, max_level,       &
                                            smth_passes, msg_fill_val)
                        else if (smth_opt == SMTHDESMTH_SPECIAL) then
                           call smth_desmth_special(field,              &
                                            start_patch_i, end_patch_i, &
                                            start_patch_j, end_patch_j, &
                                            sm1,         em1,           &
                                            sm2,         em2,           &
                                            min_level, max_level,       &
                                            smth_passes, msg_fill_val)
                       end if
  
                     else if (grid_type == 'E') then
 
                        if (trim(fieldname) == 'HGT_M' ) then
                           topo_flag_val=1.0
                           mass_flag=1.0
                        else if (trim(fieldname) == 'HGT_V') then
                           topo_flag_val=1.0
                           mass_flag=0.0
                        else
                           topo_flag_val=0.0
                        end if
  
                        if (smth_opt == ONETWOONE) then
                           call one_two_one_egrid(field,                &
                                            start_patch_i, end_patch_i, &
                                            start_patch_j, end_patch_j, &
                                            sm1,         em1,           &
                                            sm2,         em2,           &
                                            min_level, max_level,       &
                                            smth_passes, topo_flag_val, mass_flag)
                        else if (smth_opt == SMTHDESMTH) then
                           call smth_desmth_egrid(field,                &
                                            start_patch_i, end_patch_i, &
                                            start_patch_j, end_patch_j, &
                                            sm1,         em1,           &
                                            sm2,         em2,           &
                                            min_level, max_level,       &
                                            smth_passes, topo_flag_val, mass_flag)
                        else if (smth_opt == SMTHDESMTH_SPECIAL) then
                           call mprintf(.true.,WARN,'smth-desmth_special is not currently implemented for NMM. '// &
                                                    'No smoothing will be done.')
                        end if
  
                     end if

                  end if

                  call write_field(sm1, em1, sm2, em2, &
                                   min_level, max_level, trim(fieldname), datestr, real_array=field)
        
                  ! Do we calculate directional derivatives from this field?
                  call get_dfdx_name(fieldname, gradname, istatus)
                  if (istatus == 0) then
                     allocate(slp_field(sm1:em1,sm2:em2,min_level:max_level))

                     field_count = field_count + 1
                     call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
                                  i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=gradname)

                     if (grid_type == 'C') then
                        call calc_dfdx(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, &
                                       which_domain, mapfac_ptr_x, subgrid_var)
                     else if (grid_type == 'E') then
                        call calc_dfdx(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, &
                                       which_domain, subgrid_var=subgrid_var)
                     end if
                     call write_field(sm1, em1, sm2, em2, &
                                      min_level, max_level, trim(gradname), datestr, real_array=slp_field)
                     deallocate(slp_field)
                  end if
                  call get_dfdy_name(fieldname, gradname, istatus)
                  if (istatus == 0) then
                     allocate(slp_field(sm1:em1,sm2:em2,min_level:max_level))

                     field_count = field_count + 1
                     call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
                                  i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=gradname)

                     if (grid_type == 'C') then
                        call calc_dfdy(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, &
                                       which_domain, mapfac_ptr_y, subgrid_var)
                     else if (grid_type == 'E') then
                        call calc_dfdy(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, &
                                       which_domain, subgrid_var=subgrid_var)
                     end if
                     call write_field(sm1, em1, sm2, em2, &
                                      min_level, max_level, trim(gradname), datestr, real_array=slp_field)
                     deallocate(slp_field)
                  end if
        
                  deallocate(field)
      
               ! Destination field type is CATEGORICAL
               else
                  call get_max_categories(fieldname, min_category, max_category, istatus)
                  allocate(field(sm1:em1, sm2:em2, min_category:max_category))

                  ! Do we calculate a dominant category for this field?
                  call get_domcategory_name(fieldname, domname, only_save_dominant, idomcatstatus)
        
                  if (.not. only_save_dominant) then
                     field_count = field_count + 1
                     call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
                                  i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=fieldname)
                  else
                     field_count = field_count + 1
                     call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
                                  i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
                  end if

                  if ( subgrid_var ) then
                     call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
                                     sm1, em1, sm2, em2, min_category, max_category, &
                                     processed_domain, 1, sr_x=sr_x, sr_y=sr_y)
                  else
                     call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
                                     sm1, em1, sm2, em2, min_category, max_category, &
                                     processed_domain, 1, landmask=landmask)
                  end if
       
                  ! Find fractions
                  do i=sm1, em1
                     do j=sm2, em2
                        sum = 0.
                        do k=min_category,max_category
                           sum = sum + field(i,j,k)
                        end do
                        if (sum > 0.0) then
                           do k=min_category,max_category
                              field(i,j,k) = field(i,j,k) / sum
                           end do
                        else
                           do k=min_category,max_category
                              field(i,j,k) = msg_fill_val
                           end do
                        end if
                     end do
                  end do
        
                  ! If user wants to halt when a missing value is found in output field, check now
                  if (halt_on_missing) then
                     do i=sm1, em1
                        do j=sm2, em2
                           ! Only need to examine k=1
                           if (field(i,j,1) == msg_fill_val) then
                              call mprintf(.true.,ERROR,' Missing value encountered in output field. Quitting.')
                           end if
                        end do
                     end do
                  end if
        
                  ! If we should only save the dominant category, then no need to write out fractional field
                  if (.not.only_save_dominant .or. (idomcatstatus /= 0)) then
        
                     ! Finally, we may be asked to smooth the fractional field
                     call get_smooth_option(fieldname, smth_opt, smth_passes, istatus)
                     if (istatus == 0) then
                        if (grid_type == 'C') then
                           if (smth_opt == ONETWOONE) then
                              call one_two_one(field,                    &
                                             start_patch_i, end_patch_i, &
                                             start_patch_j, end_patch_j, &
                                             sm1,         em1,           &
                                             sm2,         em2,           &
                                             min_category, max_category, &
                                             smth_passes, msg_fill_val)
                           else if (smth_opt == SMTHDESMTH) then
                              call smth_desmth(field,                    &
                                             start_patch_i, end_patch_i, &
                                             start_patch_j, end_patch_j, &
                                             sm1,         em1,           &
                                             sm2,         em2,           &
                                             min_category, max_category, &
                                             smth_passes, msg_fill_val)
                           else if (smth_opt == SMTHDESMTH_SPECIAL) then
                              call smth_desmth_special(field,            &
                                             start_patch_i, end_patch_i, &
                                             start_patch_j, end_patch_j, &
                                             sm1,         em1,           &
                                             sm2,         em2,           &
                                             min_category, max_category, &
                                             smth_passes, msg_fill_val)
                           end if
                        else if (grid_type == 'E') then
                           if (smth_opt == ONETWOONE) then
                              call one_two_one_egrid(field,              &
                                             start_patch_i, end_patch_i, &
                                             start_patch_j, end_patch_j, &
                                             sm1,         em1,           &
                                             sm2,         em2,           &
                                             min_category, max_category, &
                                             smth_passes, msg_fill_val, 1.0)
                           else if (smth_opt == SMTHDESMTH) then
                              call smth_desmth_egrid(field,              &
                                             start_patch_i, end_patch_i, &
                                             start_patch_j, end_patch_j, &
                                             sm1,         em1,           &
                                             sm2,         em2,           &
                                             min_category, max_category, &
                                             smth_passes, msg_fill_val, 1.0)
                           else if (smth_opt == SMTHDESMTH_SPECIAL) then
                              call mprintf(.true.,WARN,'smth-desmth_special is not currently implemented for NMM. '// &
                                                       'No smoothing will be done.')
                           end if
                        end if
                     end if
         
                     call write_field(sm1, em1, sm2, em2, &
                                      min_category, max_category, trim(fieldname), datestr, real_array=field)
                  end if
        
                  if (idomcatstatus == 0) then
                     call mprintf(.true.,STDOUT,'  Processing %s', s1=trim(domname))
                     allocate(dominant_field(sm1:em1, sm2:em2))
 
                     if (.not. only_save_dominant) then
                        field_count = field_count + 1
                        call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
                                     i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
                     end if

                     do i=sm1, em1
                        do j=sm2, em2
                           dominant = 0.
                           dominant_field(i,j) = real(min_category-1)
                           do k=min_category,max_category
                              if (field(i,j,k) > dominant .and. field(i,j,k) /= msg_fill_val) then 
                                 dominant_field(i,j) = real(k)
                                 dominant = field(i,j,k)
             !                  else
             !                    dominant_field(i,j) = nint(msg_fill_val)
! Maybe we should put an else clause here to set the category equal to the missing fill value?
! BUG: The problem here seems to be that, when we set a fraction equal to the missing fill value
!   above, if the last fractional index we process here has been filled, we think that the dominant
!   category should be set to the missing fill value. Perhaps we could do some check to only
!   assign the msg_fill_val if no other valid category has been assigned? But this may still not
!   work if the missing fill value is something like 0.5. Somehow use bitarrays, perhaps, to remember
!   which points are missing and which just happen to have the missing fill value?
                              end if
                           end do
                           if (dominant_field(i,j) == real(min_category-1)) dominant_field(i,j) = msg_fill_val
                        end do
                     end do
                     call write_field(sm1, em1, sm2, em2, 1, 1, &
                                      trim(domname), datestr, dominant_field)
                     deallocate(dominant_field)
                  end if
       
                  deallocate(field)

               end if
      
            end if
         end if
    
      end do

      ! Close output
      call output_close()
    
      call hash_destroy(processed_fieldnames)
    
      ! Free up memory
      if (associated(xlat_array)) deallocate(xlat_array)
      if (associated(xlon_array)) deallocate(xlon_array)
      if (grid_type == 'C') then
         if (associated(xlat_array_u)) deallocate(xlat_array_u)
         if (associated(xlon_array_u)) deallocate(xlon_array_u)
         if (associated(mapfac_array_u_x)) deallocate(mapfac_array_u_x)
         if (associated(mapfac_array_u_y)) deallocate(mapfac_array_u_y)
      end if
      if (associated(xlat_array_v)) deallocate(xlat_array_v)
      if (associated(xlon_array_v)) deallocate(xlon_array_v)
      if (associated(mapfac_array_m_x)) deallocate(mapfac_array_m_x)
      if (associated(mapfac_array_m_y)) deallocate(mapfac_array_m_y)
      if (associated(mapfac_array_v_x)) deallocate(mapfac_array_v_x)
      if (associated(mapfac_array_v_y)) deallocate(mapfac_array_v_y)
      if (associated(landmask)) deallocate(landmask)
      if (associated(xlat_array_subgrid))     deallocate(xlat_array_subgrid)
      if (associated(xlon_array_subgrid))     deallocate(xlon_array_subgrid)
      if (associated(mapfac_array_x_subgrid)) deallocate(mapfac_array_x_subgrid)
      if (associated(mapfac_array_y_subgrid)) deallocate(mapfac_array_y_subgrid)

      nullify(xlat_ptr)
      nullify(xlon_ptr)
   
   end subroutine process_tile
   
   
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: calc_field
   !
   ! Purpose: This routine fills in the "field" array with interpolated source 
   !   data. When multiple resolutions of source data are available, an appropriate
   !   resolution is chosen automatically. The specified field may either be a 
   !   continuous field or a categorical field.
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   recursive subroutine calc_field(fieldname, field, xlat_array, xlon_array, istagger, &
                                   start_i, end_i, start_j, end_j, start_k, end_k, &
                                   processed_domain, ilevel, landmask, sr_x, sr_y)
   
      use bitarray_module
      use interp_module
      use llxy_module
      use misc_definitions_module
      use proc_point_module
      use queue_module
      use source_data_module
    
      implicit none
    
      ! Arguments
      integer, intent(in) :: start_i, end_i, start_j, end_j, start_k, end_k, ilevel, istagger
      real, dimension(start_i:end_i, start_j:end_j), intent(in) :: xlat_array, xlon_array
      real, dimension(start_i:end_i, start_j:end_j, start_k:end_k), intent(inout) :: field
      real, dimension(start_i:end_i, start_j:end_j), intent(in), optional :: landmask
      integer, intent(in), optional :: sr_x, sr_y
      character (len=128), intent(in) :: fieldname
      type (bitarray), intent(inout) :: processed_domain
    
      ! Local variables
      integer :: start_src_k, end_src_k
      integer :: i, j, k, ix, iy, itype
      integer :: user_iproj, istatus
      real :: mask_val
      real :: temp
      real :: scale_factor
      real :: msg_val, msg_fill_val, threshold, src_dx, src_dy, dom_dx, dom_dy
      real :: user_stand_lon, user_truelat1, user_truelat2, user_dxkm, user_dykm, &
              user_known_x, user_known_y, user_known_lat, user_known_lon
      real, pointer, dimension(:,:,:) :: data_count
      integer, pointer, dimension(:) :: interp_type
      character (len=128) :: interp_string
      type (bitarray) :: bit_domain, level_domain
      type (queue)    :: point_queue, tile_queue
      type (q_data)   :: current_pt

      ! If this is the first trip through this routine, we need to allocate the bit array that
      !  will persist through all recursive calls, tracking which grid points have been assigned
      !  a value. 
      if (ilevel == 1) call bitarray_create(processed_domain, end_i-start_i+1, end_j-start_j+1)
    
      ! Find out the projection of the data for this "priority level" (given by ilevel)
      call get_data_projection(fieldname, user_iproj, user_stand_lon, user_truelat1, user_truelat2, &
                               user_dxkm, user_dykm, user_known_x, user_known_y, user_known_lat, &
                               user_known_lon, ilevel, istatus)
    
      ! A bad status indicates that that no data for priority level ilevel is available, and thus, that
      !   no further levels will be specified. We are done processing for this level.
      if (istatus /= 0) then
         if (ilevel == 1) call bitarray_destroy(processed_domain)
         return 
      end if
    
      ! A good status indicates that there is data for this priority level, so we store the projection
      !   of that data on a stack. The projection will be on the top of the stack (and hence will be 
      !   the "active" projection) once all higher-priority levels have been processed
      call push_source_projection(user_iproj, user_stand_lon, user_truelat1, user_truelat2, &
                               user_dxkm, user_dykm, user_dykm, user_dxkm, user_known_x, user_known_y, &
                               user_known_lat, user_known_lon)
    
      ! Before proceeding with processing for this level, though, process for the next highest priority level
      !   of source data
      call calc_field(fieldname, field, xlat_array, xlon_array, istagger, start_i, end_i, &
                      start_j, end_j, start_k, end_k, processed_domain, ilevel+1, landmask, sr_x, sr_y)
    
      ! At this point, all levels of source data with higher priority have been processed, and we can assign
      !   values to all grid points that have not already been given values from higher-priority data
    
      ! Initialize point processing module
      call proc_point_init()
    
      ! Initialize queues
      call q_init(point_queue)
      call q_init(tile_queue)
    
      ! Determine whether we will be processing categorical data or continuous data
      itype = iget_source_fieldtype(fieldname, ilevel, istatus)
      call get_interp_option(fieldname, ilevel, interp_string, istatus)
      interp_type => interp_array_from_string(interp_string)

      ! Also, check whether we will be using the cell averaging interpolator for continuous fields
      if (index(interp_string,'average_gcell') /= 0 .and. itype == CONTINUOUS) then
         call get_gcell_threshold(interp_string, threshold, istatus)
         if (istatus == 0) then
            call get_source_resolution(fieldname, ilevel, src_dx, src_dy, istatus)
            if (istatus == 0) then
               call get_domain_resolution(dom_dx, dom_dy)
               if (gridtype == 'C') then
                  if (threshold*max(src_dx,src_dy)*111. <= max(dom_dx,dom_dy)/1000.) then
                     itype = SP_CONTINUOUS
                     allocate(data_count(start_i:end_i,start_j:end_j,start_k:end_k))
                     data_count = 0.
                  end if
               else if (gridtype == 'E') then
                  if (max(src_dx,src_dy) >= threshold*max(dom_dx,dom_dy)) then
                     itype = SP_CONTINUOUS
                     allocate(data_count(start_i:end_i,start_j:end_j,start_k:end_k))
                     data_count = 0.
                  end if
               end if
            end if
         end if
      end if

      call get_missing_value(fieldname, ilevel, msg_val, istatus)
      if (istatus /= 0) msg_val = NAN
      call get_missing_fill_value(fieldname, msg_fill_val, istatus)
      if (istatus /= 0) msg_fill_val = NAN
    
      call get_masked_value(fieldname, ilevel, mask_val, istatus)
      if (istatus /= 0) mask_val = -1.
    
      if (itype == CONTINUOUS .or. itype == SP_CONTINUOUS) then
         call get_source_levels(fieldname, ilevel, start_src_k, end_src_k, istatus)
         if (istatus /= 0) return
      end if
    
      ! Initialize bitarray used to track which points have been visited and assigned values while 
      !   processing *this* priority level of data
      call bitarray_create(bit_domain, end_i-start_i+1, end_j-start_j+1)
      call bitarray_create(level_domain, end_i-start_i+1, end_j-start_j+1)
    
      ! Begin by placing a point in the tile_queue
      current_pt%lat = xlat_array(start_i,start_j)
      current_pt%lon = xlon_array(start_i,start_j)
      current_pt%x = start_i
      current_pt%y = start_j
      call q_insert(tile_queue, current_pt)
    
      ! While there are still grid points in tiles that have not yet been processed
      do while (q_isdata(tile_queue))
    
         ! Take a point from the outer queue and place it in the point_queue for processing 
         current_pt = q_remove(tile_queue)
     
         ! If this level of data is categorical (i.e., is given as an array of category indices),
         !   then first try to process the entire tile in one call to accum_categorical. Any grid
         !   points that are not given values by accum_categorical and that lie within the current
         !   tile of source data are individually assigned values in the inner loop
         if (itype == CATEGORICAL) then
     
            ! Have we already visited this point? If so, this tile has already been processed by 
            !   accum_categorical. 
            if (.not. bitarray_test(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)) then
               call q_insert(point_queue, current_pt) 
               if (.not. have_processed_tile(current_pt%lat, current_pt%lon, fieldname, ilevel)) then
                  call accum_categorical(current_pt%lat, current_pt%lon, istagger, field, &
                                         start_i, end_i, start_j, end_j, start_k, end_k, &
                                         fieldname, processed_domain, level_domain, &
                                         ilevel, msg_val, mask_val, sr_x, sr_y)
! BUG: Where do we mask out those points that are on land/water when masked=land/water is set? 
               end if
               call bitarray_set(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)
            end if
     
         else if (itype == SP_CONTINUOUS) then

            ! Have we already visited this point? If so, this tile has already been processed by 
            !   accum_continuous. 
            if (.not. bitarray_test(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)) then
               call q_insert(point_queue, current_pt) 
               if (.not. have_processed_tile(current_pt%lat, current_pt%lon, fieldname, ilevel)) then
                  call accum_continuous(current_pt%lat, current_pt%lon, istagger, field, data_count, &
                                         start_i, end_i, start_j, end_j, start_k, end_k, &
                                         fieldname, processed_domain, level_domain, &
                                         ilevel, msg_val, mask_val, sr_x, sr_y)
! BUG: Where do we mask out those points that are on land/water when masked=land/water is set? 
               end if
               call bitarray_set(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)
            end if

         else if (itype == CONTINUOUS) then
     
            ! Have we already visited this point? If so, the tile containing this point has already been
            !   processed in the inner loop.
            if (.not. bitarray_test(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)) then
               call q_insert(point_queue, current_pt) 
               call bitarray_set(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)
            end if
     
         end if
     
         ! This inner loop, where all grid points contained in the current source tile are processed
         do while (q_isdata(point_queue))
            current_pt = q_remove(point_queue)
            ix = current_pt%x
            iy = current_pt%y
      
            ! Process the current point
            if (itype == CONTINUOUS .or. itype == SP_CONTINUOUS) then
      
               ! Have we already assigned this point a value from this priority level?
               if (.not. bitarray_test(level_domain, ix-start_i+1, iy-start_j+1)) then
        
                  ! If the point was already assigned a value from a higher-priority level, no 
                  !   need to assign a new value
                  if (bitarray_test(processed_domain, ix-start_i+1, iy-start_j+1)) then
                     call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
         
                  ! Otherwise, need to assign values from this level of source data if we can
                  else
                     if (present(landmask) .and. (istagger == M .or. istagger == HH)) then
                        if (landmask(ix,iy) /= mask_val) then
                           do k=start_src_k,end_src_k
                              temp = get_point(current_pt%lat, current_pt%lon, k, &
                                               fieldname, ilevel, interp_type, msg_val)
                              if (temp /= msg_val) then
                                 field(ix, iy, k) = temp
                                 call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
                                 if (itype == SP_CONTINUOUS) data_count(ix, iy, k) = 1.0
                              else
                                 field(ix, iy, k) = msg_fill_val
                              end if
                           end do
                        else
                           do k=start_k,end_k
                              field(ix,iy,k) = msg_fill_val
                           end do
                        end if
                     else
                        do k=start_src_k,end_src_k
                           temp = get_point(current_pt%lat, current_pt%lon, k, &
                                            fieldname, ilevel, interp_type, msg_val)
                           if (temp /= msg_val) then
                              field(ix, iy, k) = temp
                              call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
                              if (itype == SP_CONTINUOUS) data_count(ix, iy, k) = 1.0
                           else
                              field(ix, iy, k) = msg_fill_val
                           end if
                        end do
                     end if
                  end if
               end if
      
            else if (itype == CATEGORICAL) then
      
               ! Have we already assigned this point a value from this priority level?
               if (.not.bitarray_test(level_domain, ix-start_i+1, iy-start_j+1)) then
       
                  ! If the point was already assigned a value from a higher-priority level, no 
                  !   need to assign a new value
                  if (bitarray_test(processed_domain, ix-start_i+1, iy-start_j+1)) then
                     call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
         
                  ! Otherwise, the point was apparently not given a value when accum_categorical
                  !   was called for the current tile, and we need to assign values from this 
                  !   level of source data if we can
                  else
                     if (present(landmask) .and. (istagger == M .or. istagger == HH)) then
                        if (landmask(ix,iy) /= mask_val) then
                           temp = get_point(current_pt%lat, current_pt%lon, 1, &
                                            fieldname, ilevel, interp_type, msg_val)
      
                           do k=start_k,end_k
                              field(ix,iy,k) = 0.
                           end do
      
                           if (temp /= msg_val) then
                              if (int(temp) >= start_k .and. int(temp) <= end_k) then
                                 field(ix, iy, int(temp)) = field(ix, iy, int(temp)) + 1.
                              else
                                 call mprintf(.true.,WARN,' Attempted to assign an invalid category '// &
                                              '%i to grid point (%i, %i)', i1=int(temp), i2=ix, i3=iy)
                              end if
                              call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
                           end if
   
                        else
                           do k=start_k,end_k
                              field(ix,iy,k) = 0.
                           end do
                        end if
                     else
                        temp = get_point(current_pt%lat, current_pt%lon, 1, &
                                         fieldname, ilevel, interp_type, msg_val)
      
                        do k=start_k,end_k
                           field(ix,iy,k) = 0.
                        end do
      
                        if (temp /= msg_val) then
                           if (int(temp) >= start_k .and. int(temp) <= end_k) then
                              field(ix, iy, int(temp)) = field(ix, iy, int(temp)) + 1.
                           else
                              call mprintf(.true.,WARN,' Attempted to assign an invalid category '// &
                                           '%i to grid point (%i, %i)', i1=int(temp), i2=ix, i3=iy)
                           end if
                           call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
                        end if
                     end if
                  end if
               end if
      
            end if
      
            ! Scan neighboring points, adding them to the appropriate queue based on whether they
            !   are in the current tile or not
            if (iy > start_j) then
               if (ix > start_i) then
        
                  ! Neighbor with relative position (-1,-1)
                  call process_neighbor(ix-1, iy-1, bit_domain, point_queue, tile_queue, &
                                        xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
               end if
       
               ! Neighbor with relative position (0,-1) 
               call process_neighbor(ix, iy-1, bit_domain, point_queue, tile_queue, &
                                     xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
       
               if (ix < end_i) then
       
                  ! Neighbor with relative position (+1,-1)
                  call process_neighbor(ix+1, iy-1, bit_domain, point_queue, tile_queue, &
                                        xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
               end if
            end if
      
            if (ix > start_i) then
      
               ! Neighbor with relative position (-1,0)
               call process_neighbor(ix-1, iy, bit_domain, point_queue, tile_queue, &
                                     xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
            end if
      
            if (ix < end_i) then
      
               ! Neighbor with relative position (+1,0)
               call process_neighbor(ix+1, iy, bit_domain, point_queue, tile_queue, &
                                     xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
            end if
      
            if (iy < end_j) then
               if (ix > start_i) then
       
                  ! Neighbor with relative position (-1,+1)
                  call process_neighbor(ix-1, iy+1, bit_domain, point_queue, tile_queue, &
                                        xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
               end if
       
               ! Neighbor with relative position (0,+1)
               call process_neighbor(ix, iy+1, bit_domain, point_queue, tile_queue, &
                                     xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
               if (ix < end_i) then
       
                  ! Neighbor with relative position (+1,+1)
                  call process_neighbor(ix+1, iy+1, bit_domain, point_queue, tile_queue, &
                                        xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
               end if
            end if
      
         end do
     
      end do

      if (itype == SP_CONTINUOUS) then
         itype = CONTINUOUS
         if (present(landmask) .and. (istagger == M .or. istagger == HH)) then
            do j=start_j,end_j
               do i=start_i,end_i
                  if (landmask(i,j) /= mask_val) then
                     do k=start_k,end_k
                        if (data_count(i,j,k) > 0.) then
                           field(i,j,k) = field(i,j,k) / data_count(i,j,k)
                        else
                           if (.not.bitarray_test(processed_domain, i-start_i+1, j-start_j+1)) then
                              field(i,j,k) = msg_fill_val
                           end if
                        end if
                     end do
                  else
                     if (.not.bitarray_test(processed_domain, i-start_i+1, j-start_j+1)) then
                        do k=start_k,end_k
                           field(i,j,k) = msg_fill_val
                        end do
                     end if
                  end if
               end do
            end do
         else
            do k=start_k,end_k
               do j=start_j,end_j
                  do i=start_i,end_i
                     if (data_count(i,j,k) > 0.) then
                        field(i,j,k) = field(i,j,k) / data_count(i,j,k)
                     else
                        if (.not.bitarray_test(processed_domain, i-start_i+1, j-start_j+1)) then
                           field(i,j,k) = msg_fill_val
                        end if
                     end if
                  end do
               end do
            end do
         end if
         deallocate(data_count)

      else if (itype == CATEGORICAL) then
         if (present(landmask) .and. (istagger == M .or. istagger == HH)) then
            do j=start_j,end_j
               do i=start_i,end_i
                  if (landmask(i,j) == mask_val) then
                     do k=start_k,end_k
                        field(i,j,k) = 0.
                     end do
                  end if
               end do
            end do
         end if
      end if

      deallocate(interp_type)


      ! We may need to scale this field by a constant
      call get_field_scale_factor(fieldname, ilevel, scale_factor, istatus)
      if (istatus == 0) then
         do i=start_i, end_i
            do j=start_j, end_j
               if (bitarray_test(level_domain,i-start_i+1,j-start_j+1) .and. &
                   .not. bitarray_test(processed_domain,i-start_i+1,j-start_j+1)) then
                  do k=start_k,end_k
                     if (field(i,j,k) /= msg_fill_val) then
                        field(i,j,k) = field(i,j,k) * scale_factor
                     end if
                  end do
               end if
            end do
         end do
      end if

    
      ! Now add the points that were assigned values at this priority level to the complete array
      !   of points that have been assigned values
      call bitarray_merge(processed_domain, level_domain)
    
      call bitarray_destroy(bit_domain)
      call bitarray_destroy(level_domain)
      call q_destroy(point_queue)
      call q_destroy(tile_queue)
      call proc_point_shutdown()
      ! If this is the last level of the recursion, we can also deallocate processed_domain
      if (ilevel == 1) call bitarray_destroy(processed_domain)
    
      ! Remove the projection of the current level of source data from the stack, thus "activating" 
      !   the projection of the next highest level
      call pop_source_projection()
   
   end subroutine calc_field
   
   
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: get_lat_lon_fields
   !
   ! Purpose: To calculate the latitude and longitude for every gridpoint in the
   !   tile of the model domain. The caller may specify that the grid for which 
   !   values are computed is staggered or unstaggered using the "stagger" 
   !   argument.
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   subroutine get_lat_lon_fields(xlat_arr, xlon_arr, start_mem_i, &
                                 start_mem_j, end_mem_i, end_mem_j, stagger, comp_ll, &
                                 sub_x, sub_y)
   
      use llxy_module
      use misc_definitions_module
    
      implicit none
    
      ! Arguments
      integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, &
                             end_mem_j, stagger
      real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: xlat_arr, xlon_arr
      logical, optional, intent(in) :: comp_ll
      integer, optional, intent(in) :: sub_x, sub_y

      ! Local variables
      integer :: i, j
      real :: rx, ry, dx, dy
    
      rx = 1.0
      ry = 1.0
      dx = 1.0
      dy = 1.0

      if (present(sub_x)) then 
          rx = real(sub_x)
          dx = .5 + .5/rx
      end if
      if (present(sub_y)) then
          ry = real(sub_y)
          dy = .5 + .5/ry
      end if

      do i=start_mem_i, end_mem_i
         do j=start_mem_j, end_mem_j
            call xytoll(real(i-1)/rx+dx, real(j-1)/ry+dy, &
                        xlat_arr(i,j), xlon_arr(i,j), stagger, comp_ll=comp_ll)
         end do
      end do

   end subroutine get_lat_lon_fields
   

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: get_map_factor
   !
   ! Purpose: Given the latitude field, this routine calculates map factors for 
   !   the grid points of the specified domain. For different grids (e.g., C grid, 
   !   E grid), the latitude array should provide the latitudes of the points for
   !   which map factors are to be calculated. 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   subroutine get_map_factor(xlat_arr, xlon_arr, mapfac_arr_x, mapfac_arr_y, &
                             start_mem_i, start_mem_j, end_mem_i, end_mem_j)
   
      use constants_module
      use gridinfo_module
      use misc_definitions_module
      use map_utils
    
      implicit none
    
      ! Arguments
      integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, end_mem_j
      real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in) :: xlat_arr, xlon_arr
      real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: mapfac_arr_x
      real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: mapfac_arr_y
    
      ! Local variables
      integer :: i, j
      real :: n, colat, colat0, colat1, colat2, comp_lat, comp_lon
    
      !
      ! Equations for map factor given in Principles of Meteorological Analysis,
      ! Walter J. Saucier, pp. 32-33 
      !
    
      ! Lambert conformal projection
      if (iproj_type == PROJ_LC) then
         if (truelat1 /= truelat2) then
            colat1 = rad_per_deg*(90.0 - truelat1)
            colat2 = rad_per_deg*(90.0 - truelat2)
            n = (log(sin(colat1)) - log(sin(colat2))) &
                / (log(tan(colat1/2.0)) - log(tan(colat2/2.0)))
      
            do i=start_mem_i, end_mem_i
               do j=start_mem_j, end_mem_j
                  colat = rad_per_deg*(90.0 - xlat_arr(i,j))
                  mapfac_arr_x(i,j) = sin(colat2)/sin(colat)*(tan(colat/2.0)/tan(colat2/2.0))**n
                  mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
               end do
            end do
     
         else
            colat0 = rad_per_deg*(90.0 - truelat1)
      
            do i=start_mem_i, end_mem_i
               do j=start_mem_j, end_mem_j
                  colat = rad_per_deg*(90.0 - xlat_arr(i,j))
                  mapfac_arr_x(i,j) = sin(colat0)/sin(colat)*(tan(colat/2.0)/tan(colat0/2.0))**cos(colat0)
                  mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
               end do
            end do
    
         end if
    
      ! Polar stereographic projection
      else if (iproj_type == PROJ_PS) then
    
         do i=start_mem_i, end_mem_i
            do j=start_mem_j, end_mem_j
               mapfac_arr_x(i,j) = (1.0 + sin(rad_per_deg*abs(truelat1)))/(1.0 + sin(rad_per_deg*sign(1.,truelat1)*xlat_arr(i,j)))
               mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
            end do
         end do
    
      ! Mercator projection 
      else if (iproj_type == PROJ_MERC) then
         colat0 = rad_per_deg*(90.0 - truelat1)
     
         do i=start_mem_i, end_mem_i
            do j=start_mem_j, end_mem_j
               colat = rad_per_deg*(90.0 - xlat_arr(i,j))
               mapfac_arr_x(i,j) = sin(colat0) / sin(colat) 
               mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
            end do
         end do
    
      ! Global cylindrical projection
      else if (iproj_type == PROJ_CYL) then
     
         do i=start_mem_i, end_mem_i
            do j=start_mem_j, end_mem_j
               if (abs(xlat_arr(i,j)) == 90.0) then
                  mapfac_arr_x(i,j) = 0.    ! MSF actually becomes infinite at poles, but 
                                            !   the values should never be used there; by
                                            !   setting to 0, we hope to induce a "divide
                                            !   by zero" error if they are
               else
                  mapfac_arr_x(i,j) = 1.0 / cos(xlat_arr(i,j)*rad_per_deg) 
               end if
               mapfac_arr_y(i,j) = 1.0
            end do
         end do
    
      ! Rotated global cylindrical projection
      else if (iproj_type == PROJ_CASSINI) then
     
         if (abs(pole_lat) == 90.) then
            do i=start_mem_i, end_mem_i
               do j=start_mem_j, end_mem_j
                  if (abs(xlat_arr(i,j)) >= 90.0) then
                     mapfac_arr_x(i,j) = 0.    ! MSF actually becomes infinite at poles, but 
                                               !   the values should never be used there; by
                                               !   setting to 0, we hope to induce a "divide
                                               !   by zero" error if they are
                  else
                     mapfac_arr_x(i,j) = 1.0 / cos(xlat_arr(i,j)*rad_per_deg) 
                  end if
                  mapfac_arr_y(i,j) = 1.0
               end do
            end do
         else
            do i=start_mem_i, end_mem_i
               do j=start_mem_j, end_mem_j
                  call rotate_coords(xlat_arr(i,j),xlon_arr(i,j), &
                                     comp_lat, comp_lon, &
                                     pole_lat, pole_lon, stand_lon, &
                                     -1)
                  if (abs(comp_lat) >= 90.0) then
                     mapfac_arr_x(i,j) = 0.    ! MSF actually becomes infinite at poles, but 
                                               !   the values should never be used there; by
                                               !   setting to 0, we hope to induce a "divide
                                               !   by zero" error if they are
                  else
                     mapfac_arr_x(i,j) = 1.0 / cos(comp_lat*rad_per_deg) 
                  end if
                  mapfac_arr_y(i,j) = 1.0
               end do
            end do
         end if
    
      else if (iproj_type == PROJ_ROTLL) then
    
         do i=start_mem_i, end_mem_i
            do j=start_mem_j, end_mem_j
               mapfac_arr_x(i,j) = 1.0
               mapfac_arr_y(i,j) = 1.0
            end do
         end do
    
      end if
   
   end subroutine get_map_factor
   
   
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: get_coriolis_parameters
   !
   ! Purpose: To calculate the Coriolis parameters f and e for every gridpoint in
   !   the tile of the model domain 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   subroutine get_coriolis_parameters(xlat_arr, f, e, &
                                      start_mem_i, start_mem_j, end_mem_i, end_mem_j)
     
      use constants_module
    
      implicit none
    
      ! Arguments
      integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, end_mem_j
      real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in) :: xlat_arr
      real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: f, e
    
      ! Local variables
      integer :: i, j
    
      do i=start_mem_i, end_mem_i
         do j=start_mem_j, end_mem_j
     
            f(i,j) = 2.0*OMEGA_E*sin(rad_per_deg*xlat_arr(i,j))
            e(i,j) = 2.0*OMEGA_E*cos(rad_per_deg*xlat_arr(i,j))
     
         end do
      end do
      
   end subroutine get_coriolis_parameters
   
   
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: get_rotang
   !
   ! Purpose: To calculate the sine and cosine of rotation angle. 
   !
   ! NOTES: The formulas used in this routine come from those in the 
   !   vecrot_rotlat() routine of the original WRF SI.
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   subroutine get_rotang(xlat_arr, xlon_arr, cosa, sina, &
                         start_mem_i, start_mem_j, end_mem_i, end_mem_j)
   
      use constants_module
      use gridinfo_module
    
      implicit none
    
      ! Arguments
      integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, end_mem_j
      real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in) :: xlat_arr, xlon_arr
      real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: cosa, sina
    
      ! Local variables
      integer :: i, j
      real :: alpha, d_lon

      do i=start_mem_i, end_mem_i
         do j=start_mem_j+1, end_mem_j-1
            d_lon = xlon_arr(i,j+1)-xlon_arr(i,j-1)
            if (d_lon > 180.) then
               d_lon = d_lon - 360.
            else if (d_lon < -180.) then
               d_lon = d_lon + 360.
            end if

            alpha = atan2(-cos(xlat_arr(i,j)*RAD_PER_DEG) * (d_lon*RAD_PER_DEG), &
                            ((xlat_arr(i,j+1)-xlat_arr(i,j-1))*RAD_PER_DEG))
            sina(i,j) = sin(alpha)
            cosa(i,j) = cos(alpha)
         end do
      end do

      do i=start_mem_i, end_mem_i
         d_lon = xlon_arr(i,start_mem_j+1)-xlon_arr(i,start_mem_j)
         if (d_lon > 180.) then
            d_lon = d_lon - 360.
         else if (d_lon < -180.) then
            d_lon = d_lon + 360.
         end if

         alpha = atan2(-cos(xlat_arr(i,start_mem_j)*RAD_PER_DEG) * (d_lon*RAD_PER_DEG), &
                       ((xlat_arr(i,start_mem_j+1)-xlat_arr(i,start_mem_j))*RAD_PER_DEG))
         sina(i,start_mem_j) = sin(alpha)
         cosa(i,start_mem_j) = cos(alpha)
      end do

      do i=start_mem_i, end_mem_i
         d_lon = xlon_arr(i,end_mem_j)-xlon_arr(i,end_mem_j-1)
         if (d_lon > 180.) then
            d_lon = d_lon - 360.
         else if (d_lon < -180.) then
            d_lon = d_lon + 360.
         end if

         alpha = atan2(-cos(xlat_arr(i,end_mem_j)*RAD_PER_DEG) * (d_lon*RAD_PER_DEG), &
                       ((xlat_arr(i,end_mem_j)-xlat_arr(i,end_mem_j-1))*RAD_PER_DEG))
         sina(i,end_mem_j) = sin(alpha)
         cosa(i,end_mem_j) = cos(alpha)
      end do
    
   end subroutine get_rotang
   
   
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: process_neighbor
   !
   ! Purpose: This routine, give the x/y location of a point, determines whether
   !   the point has already been processed, and if not, which processing queue
   !   the point should be placed in.
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   subroutine process_neighbor(ix, iy, bit_domain, point_queue, tile_queue, &
                               xlat_array, xlon_array, &
                               start_i, end_i, start_j, end_j, ilevel)
   
      use bitarray_module
      use misc_definitions_module
      use proc_point_module
      use queue_module
    
      implicit none
    
      ! Arguments
      integer, intent(in) :: ix, iy, start_i, end_i, start_j, end_j, ilevel
      real, dimension(start_i:end_i, start_j:end_j), intent(in) :: xlat_array, xlon_array
      type (bitarray), intent(inout) :: bit_domain
      type (queue), intent(inout) :: point_queue, tile_queue
    
      ! Local variables
      type (q_data) :: process_pt
      logical :: is_in_tile
    
      ! If the point has already been visited, no need to do anything more.
      if (.not. bitarray_test(bit_domain, ix-start_i+1, iy-start_j+1)) then 
    
         ! Create a queue item for the current point
         process_pt%lat = xlat_array(ix,iy)            
         process_pt%lon = xlon_array(ix,iy)            
         process_pt%x = ix
         process_pt%y = iy
     
         is_in_tile = is_point_in_tile(process_pt%lat, process_pt%lon, ilevel)
     
         ! If the point is in the current tile, add it to the list of points
         !   to be processed in the inner loop
         if (is_in_tile) then
            call q_insert(point_queue, process_pt)
            call bitarray_set(bit_domain, ix-start_i+1, iy-start_j+1) 
     
         ! Otherwise, we will process this point later. Add it to the list for 
         !   the outer loop.
         else
            call q_insert(tile_queue, process_pt)
         end if
    
      end if
   
   end subroutine process_neighbor
   
   
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: calc_dfdy
   !
   ! Purpose: This routine calculates df/dy for the field in src_arr, and places
   !   the result in dst_array.
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   subroutine calc_dfdy(src_arr, dst_arr, start_mem_i, start_mem_j, start_mem_k, &
                        end_mem_i, end_mem_j, end_mem_k, idom, mapfac, subgrid_var)
   
      ! Modules
      use gridinfo_module
    
      implicit none
    
      ! Arguments
      integer, intent(in) :: start_mem_i, start_mem_j, start_mem_k, end_mem_i, end_mem_j, end_mem_k
      real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j,start_mem_k:end_mem_k), intent(in) :: src_arr
      real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j,start_mem_k:end_mem_k), intent(out) :: dst_arr
      integer,intent(in) :: idom
      real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in), optional :: mapfac
      logical, intent(in), optional :: subgrid_var
    
      ! Local variables
      integer :: i, j, k
      real :: l_sr_y

      ! calculate the refinement ratio from the top-level domain
      l_sr_y=1.
      if (present(subgrid_var)) then
        if (subgrid_var) l_sr_y=subgrid_ratio_y(idom)
      end if
      l_sr_y=l_sr_y*moad_grid_ratio(idom)

      if (present(mapfac)) then
         do k=start_mem_k,end_mem_k
            do i=start_mem_i, end_mem_i
               do j=start_mem_j+1, end_mem_j-1
                  dst_arr(i,j,k) = l_sr_y*(src_arr(i,j+1,k) - src_arr(i,j-1,k))/(2.*dykm*mapfac(i,j))
               end do
            end do
     
            do i=start_mem_i, end_mem_i
               dst_arr(i,start_mem_j,k) = l_sr_y*(src_arr(i,start_mem_j+1,k) - src_arr(i,start_mem_j,k))/(dykm*mapfac(i,j))
            end do
     
            do i=start_mem_i, end_mem_i
               dst_arr(i,end_mem_j,k) = l_sr_y*(src_arr(i,end_mem_j,k) - src_arr(i,end_mem_j-1,k))/(dykm*mapfac(i,j))
            end do
         end do
      else
         do k=start_mem_k,end_mem_k
            do i=start_mem_i, end_mem_i
               do j=start_mem_j+1, end_mem_j-1
                  dst_arr(i,j,k) = l_sr_y*(src_arr(i,j+1,k) - src_arr(i,j-1,k))/(2.*dykm)
               end do
            end do
     
            do i=start_mem_i, end_mem_i
               dst_arr(i,start_mem_j,k) = l_sr_y*(src_arr(i,start_mem_j+1,k) - src_arr(i,start_mem_j,k))/(dykm)
            end do
     
            do i=start_mem_i, end_mem_i
               dst_arr(i,end_mem_j,k) = l_sr_y*(src_arr(i,end_mem_j,k) - src_arr(i,end_mem_j-1,k))/(dykm)
            end do
         end do
      end if
   
   end subroutine calc_dfdy
   
   
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: calc_dfdx
   !
   ! Purpose: This routine calculates df/dx for the field in src_arr, and places
   !   the result in dst_array.
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   subroutine calc_dfdx(src_arr, dst_arr, start_mem_i, start_mem_j, &
                        start_mem_k, end_mem_i, end_mem_j, end_mem_k, idom, mapfac, subgrid_var)
   
      ! Modules
      use gridinfo_module
    
      implicit none
    
      ! Arguments
      integer, intent(in) :: start_mem_i, start_mem_j, start_mem_k, end_mem_i, end_mem_j, end_mem_k
      real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j, start_mem_k:end_mem_k), intent(in) :: src_arr
      real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j, start_mem_k:end_mem_k), intent(out) :: dst_arr
      integer, intent(in) :: idom
      real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in), optional :: mapfac
      logical, intent(in), optional :: subgrid_var
    
      ! Local variables
      integer :: i, j, k
      real :: l_sr_x

      ! calculate the refinement ratio from the top-level domain
      l_sr_x=1.
      if (present(subgrid_var)) then 
        if (subgrid_var) l_sr_x=subgrid_ratio_x(idom)
      end if
      l_sr_x=l_sr_x*moad_grid_ratio(idom)

      if (present(mapfac)) then
         do k=start_mem_k, end_mem_k
            do i=start_mem_i+1, end_mem_i-1
               do j=start_mem_j, end_mem_j
                  dst_arr(i,j,k) = l_sr_x*(src_arr(i+1,j,k) - src_arr(i-1,j,k))/(2.*dxkm*mapfac(i,j))
               end do
            end do
     
            do j=start_mem_j, end_mem_j
               dst_arr(start_mem_i,j,k) = l_sr_x*(src_arr(start_mem_i+1,j,k) - src_arr(start_mem_i,j,k))/(dxkm*mapfac(i,j))
            end do
     
            do j=start_mem_j, end_mem_j
               dst_arr(end_mem_i,j,k) = l_sr_x*(src_arr(end_mem_i,j,k) - src_arr(end_mem_i-1,j,k))/(dxkm*mapfac(i,j))
            end do
         end do
      else
         do k=start_mem_k, end_mem_k
            do i=start_mem_i+1, end_mem_i-1
               do j=start_mem_j, end_mem_j
                  dst_arr(i,j,k) = l_sr_x*(src_arr(i+1,j,k) - src_arr(i-1,j,k))/(2.*dxkm)
               end do
            end do
     
            do j=start_mem_j, end_mem_j
               dst_arr(start_mem_i,j,k) = l_sr_x*(src_arr(start_mem_i+1,j,k) - src_arr(start_mem_i,j,k))/(dxkm)
            end do
     
            do j=start_mem_j, end_mem_j
               dst_arr(end_mem_i,j,k) = l_sr_x*(src_arr(end_mem_i,j,k) - src_arr(end_mem_i-1,j,k))/(dxkm)
            end do
         end do
      end if
   
   end subroutine calc_dfdx

end module process_tile_module