metgrid/src/process_domain_module.F

   1 module process_domain_module
   2
   3    use mpas_mesh
   4    use target_mesh
   5    use remapper
   6
   7    type (mpas_mesh_type), save :: mpas_source_mesh
   8    type (target_mesh_type), save :: wrf_target_mesh_m, wrf_target_mesh_u, wrf_target_mesh_v
   9    type (remap_info_type), save :: remap_info_m, remap_info_u, remap_info_v
  10    real, dimension(:,:), allocatable, target :: xlat_rad, xlon_rad, xlat_u_rad, xlon_u_rad, xlat_v_rad, xlon_v_rad
  11
  12    contains
  13
  14    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  15    ! Name: process_domain
  16    !
  17    ! Purpose:
  18    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  19    subroutine process_domain(n, extra_row, extra_col)
  20
  21       use date_pack
  22       use gridinfo_module
  23       use interp_option_module
  24       use misc_definitions_module
  25       use module_debug
  26       use storage_module
  27
  28       implicit none
  29
  30       ! Arguments
  31       integer, intent(in) :: n
  32       logical, intent(in) :: extra_row, extra_col
  33
  34       ! Local variables
  35       integer :: istatus
  36       integer :: i, t, dyn_opt, &
  37                  we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
  38                  we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
  39                  sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
  40                  we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
  41                  sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
  42                  idiff, n_times, &
  43                  west_east_dim, south_north_dim, bottom_top_dim, map_proj, &
  44                  is_water, is_lake, is_ice, is_urban, i_soilwater, &
  45                  grid_id, parent_id, i_parent_start, j_parent_start, &
  46                  i_parent_end, j_parent_end, parent_grid_ratio, sub_x, sub_y, num_land_cat, process_bdy_width
  47       real :: cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
  48               dom_dx, dom_dy, pole_lat, pole_lon
  49       real, dimension(16) :: corner_lats, corner_lons
  50       real, pointer, dimension(:,:) :: landmask
  51       real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
  52       logical, allocatable, dimension(:) :: got_this_field, got_const_field
  53       character (len=19) :: valid_date, temp_date
  54       character (len=128) :: title, mminlu
  55       character (len=128), allocatable, dimension(:) :: output_flags, td_output_flags
  56       character (len=128), dimension(:), pointer :: geogrid_flags
  57
  58       ! CWH Initialize local pointer variables
  59       nullify(landmask)
  60       nullify(xlat)
  61       nullify(xlon)
  62       nullify(xlat_u)
  63       nullify(xlon_u)
  64       nullify(xlat_v)
  65       nullify(xlon_v)
  66       nullify(geogrid_flags)
  67
  68       ! Compute number of times that we will process
  69       call geth_idts(end_date(n), start_date(n), idiff)
  70       call mprintf((idiff < 0),ERROR,'Ending date is earlier than starting date in namelist for domain %i.', i1=n)
  71
  72       n_times = idiff / interval_seconds
  73
  74       ! Check that the interval evenly divides the range of times to process
  75       call mprintf((mod(idiff, interval_seconds) /= 0),WARN, &
  76                    'In namelist, interval_seconds does not evenly divide '// &
  77                    '(end_date - start_date) for domain %i. Only %i time periods '// &
  78                    'will be processed.', i1=n, i2=n_times)
  79
  80       ! Initialize the storage module
  81       call mprintf(.true.,LOGFILE,'Initializing storage module')
  82       call storage_init()
  83
  84       !
  85       ! Do time-independent processing
  86       !
  87       call get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, map_proj, &
  88                     we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
  89                     we_patch_s,      we_patch_e, &
  90                     we_patch_stag_s, we_patch_stag_e, &
  91                     sn_patch_s,      sn_patch_e, &
  92                     sn_patch_stag_s, sn_patch_stag_e, &
  93                     we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
  94                     sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
  95                     mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
  96                     grid_id, parent_id, &
  97                     i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
  98                     parent_grid_ratio, sub_x, sub_y, &
  99                     cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
 100                     pole_lat, pole_lon, dom_dx, dom_dy, landmask, xlat, xlon, xlat_u, xlon_u, &
 101                     xlat_v, xlon_v, corner_lats, corner_lons, title, geogrid_flags)
 102
 103
 104       allocate(output_flags(num_entries))
 105       allocate(got_const_field(num_entries))
 106
 107       do i=1,num_entries
 108          output_flags(i)    = ' '
 109          got_const_field(i) = .false.
 110       end do
 111
 112       ! This call is to process the constant met fields (SST or SEAICE, for example)
 113       ! That we process constant fields is indicated by the first argument
 114       call process_single_met_time(.true., temp_date, n, extra_row, extra_col, xlat, xlon, &
 115                           xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
 116                           title, dyn_opt, &
 117                           west_east_dim, south_north_dim, &
 118                           we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
 119                           we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
 120                           sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
 121                           we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
 122                           sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
 123                           got_const_field, &
 124                           map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
 125                           grid_id, parent_id, i_parent_start, &
 126                           j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
 127                           cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
 128                           pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
 129                           corner_lats, corner_lons, output_flags, geogrid_flags, 0)
 130
 131       !
 132       ! Begin time-dependent processing
 133       !
 134
 135       allocate(td_output_flags(num_entries))
 136       allocate(got_this_field (num_entries))
 137
 138       ! Loop over all times to be processed for this domain
 139       do t=0,n_times
 140
 141          call geth_newdate(valid_date, trim(start_date(n)), t*interval_seconds)
 142          temp_date = ' '
 143
 144          if (mod(interval_seconds,3600) == 0) then
 145             write(temp_date,'(a13)') valid_date(1:10)//'_'//valid_date(12:13)
 146          else if (mod(interval_seconds,60) == 0) then
 147             write(temp_date,'(a16)') valid_date(1:10)//'_'//valid_date(12:16)
 148          else
 149             write(temp_date,'(a19)') valid_date(1:10)//'_'//valid_date(12:19)
 150          end if
 151
 152          call mprintf(.true.,STDOUT, ' Processing %s', s1=trim(temp_date))
 153          call mprintf(.true.,LOGFILE, 'Preparing to process output time %s', s1=temp_date)
 154
 155          do i=1,num_entries
 156             td_output_flags(i) = output_flags(i)
 157             got_this_field(i)  = got_const_field(i)
 158          end do
 159
 160          if (t > 0) then
 161             process_bdy_width = process_only_bdy
 162          else
 163             process_bdy_width = 0
 164          end if
 165
 166          call process_single_met_time(.false., temp_date, n, extra_row, extra_col, xlat, xlon, &
 167                              xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
 168                              title, dyn_opt, &
 169                              west_east_dim, south_north_dim, &
 170                              we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
 171                              we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
 172                              sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
 173                              we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
 174                              sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
 175                              got_this_field, &
 176                              map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
 177                              grid_id, parent_id, i_parent_start, &
 178                              j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
 179                              cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
 180                              pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
 181                              corner_lats, corner_lons, td_output_flags, geogrid_flags, process_bdy_width)
 182
 183       end do  ! Loop over n_times
 184
 185
 186       deallocate(td_output_flags)
 187       deallocate(got_this_field)
 188
 189       deallocate(output_flags)
 190       deallocate(got_const_field)
 191
 192       if (associated(geogrid_flags)) deallocate(geogrid_flags)
 193
 194       call storage_delete_all()
 195
 196       istatus = mpas_mesh_free(mpas_source_mesh)
 197
 198       if (allocated(xlat_rad)) deallocate(xlat_rad)
 199       if (allocated(xlon_rad)) deallocate(xlon_rad)
 200       if (allocated(xlat_u_rad)) deallocate(xlat_u_rad)
 201       if (allocated(xlon_u_rad)) deallocate(xlon_u_rad)
 202       if (allocated(xlat_v_rad)) deallocate(xlat_v_rad)
 203       if (allocated(xlon_v_rad)) deallocate(xlon_v_rad)
 204       istatus = target_mesh_free(wrf_target_mesh_m)
 205       istatus = target_mesh_free(wrf_target_mesh_u)
 206       istatus = target_mesh_free(wrf_target_mesh_v)
 207
 208       istatus = remap_info_free(remap_info_m)
 209       istatus = remap_info_free(remap_info_u)
 210       istatus = remap_info_free(remap_info_v)
 211
 212    end subroutine process_domain
 213
 214
 215    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 216    ! Name: get_static_fields
 217    !
 218    ! Purpose:
 219    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 220    subroutine get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, &
 221                     map_proj,                                                               &
 222                     we_dom_s,   we_dom_e,   sn_dom_s,        sn_dom_e,                      &
 223                     we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e,               &
 224                     sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e,               &
 225                     we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e,                       &
 226                     sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e,                       &
 227                     mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
 228                     grid_id, parent_id,                                                     &
 229                     i_parent_start, j_parent_start, i_parent_end, j_parent_end,             &
 230                     parent_grid_ratio, sub_x, sub_y,                                        &
 231                     cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2,          &
 232                     pole_lat, pole_lon, dom_dx, dom_dy, landmask, xlat, xlon, xlat_u, xlon_u, &
 233                     xlat_v, xlon_v, corner_lats, corner_lons, title, geogrid_flags)
 234
 235       use gridinfo_module
 236       use input_module
 237       use llxy_module
 238       use parallel_module
 239       use storage_module
 240       use module_debug
 241       use list_module
 242
 243       implicit none
 244
 245       ! Arguments
 246       integer, intent(in) :: n
 247       integer, intent(inout) :: dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, &
 248                                 map_proj, &
 249                                 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
 250                                 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
 251                                 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
 252                                 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
 253                                 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
 254                                 is_water, is_lake, is_ice, is_urban, i_soilwater, grid_id, parent_id, &
 255                                 i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
 256                                 parent_grid_ratio, sub_x, sub_y, num_land_cat
 257       real, pointer, dimension(:,:) :: landmask
 258       real, intent(inout) :: cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
 259                              dom_dx, dom_dy, pole_lat, pole_lon
 260       real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
 261       real, dimension(16), intent(out) :: corner_lats, corner_lons
 262       character (len=128), intent(inout) :: title, mminlu
 263       character (len=128), dimension(:), pointer :: geogrid_flags
 264
 265       ! Local variables
 266       integer :: istatus, i, j, k, sp1, ep1, sp2, ep2, sp3, ep3, &
 267                  lh_mult, rh_mult, bh_mult, th_mult, subx, suby
 268       integer :: we_mem_subgrid_s, we_mem_subgrid_e, &
 269                  sn_mem_subgrid_s, sn_mem_subgrid_e
 270       integer :: we_patch_subgrid_s, we_patch_subgrid_e, &
 271                  sn_patch_subgrid_s, sn_patch_subgrid_e
 272       real, pointer, dimension(:,:,:) :: real_array
 273       character (len=3) :: memorder
 274       character (len=128) :: grid_type, datestr, cname, stagger, cunits, cdesc, ctemp
 275       character (len=128), dimension(3) :: dimnames
 276       type (fg_input) :: field
 277       type (list) :: static_list, flag_list
 278       type (list_item), dimension(:), pointer :: static_list_array
 279
 280       call list_init(static_list)
 281
 282       ! CWH Initialize local pointer variables
 283       nullify(real_array)
 284
 285       ! Initialize the input module to read static input data for this domain
 286       call mprintf(.true.,LOGFILE,'Opening static input file.')
 287       call input_init(n, istatus)
 288       call mprintf((istatus /= 0),ERROR, 'input_init(): Error opening input for domain %i.', i1=n)
 289
 290       ! Read global attributes from the static data input file
 291       call mprintf(.true.,LOGFILE,'Reading static global attributes.')
 292       call read_global_attrs(title, datestr, grid_type, dyn_opt, west_east_dim,          &
 293                              south_north_dim, bottom_top_dim,                            &
 294                              we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e,   &
 295                              sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e,   &
 296                              map_proj, mminlu, num_land_cat,                             &
 297                              is_water, is_lake, is_ice, is_urban, i_soilwater,           &
 298                              grid_id, parent_id, i_parent_start,                         &
 299                              j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
 300                              cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1,        &
 301                              truelat2, pole_lat, pole_lon, parent_grid_ratio,            &
 302                              corner_lats, corner_lons, sub_x, sub_y)
 303
 304       we_dom_s = 1
 305       sn_dom_s = 1
 306       if (grid_type(1:1) == 'C') then
 307          we_dom_e = west_east_dim   - 1
 308          sn_dom_e = south_north_dim - 1
 309       else if (grid_type(1:1) == 'E') then
 310          we_dom_e = west_east_dim
 311          sn_dom_e = south_north_dim
 312       end if
 313
 314       ! Given the full dimensions of this domain, find out the range of indices
 315       !   that will be worked on by this processor. This information is given by
 316       !   my_minx, my_miny, my_maxx, my_maxy
 317       call parallel_get_tile_dims(west_east_dim, south_north_dim)
 318
 319       ! Must figure out patch dimensions from info in parallel module
 320       if (nprocs > 1 .and. .not. do_tiled_input) then
 321
 322          we_patch_s      = my_minx
 323          we_patch_stag_s = my_minx
 324          we_patch_e      = my_maxx - 1
 325          sn_patch_s      = my_miny
 326          sn_patch_stag_s = my_miny
 327          sn_patch_e      = my_maxy - 1
 328
 329          if (gridtype == 'C') then
 330             if (my_x /= nproc_x - 1) then
 331                we_patch_e = we_patch_e + 1
 332                we_patch_stag_e = we_patch_e
 333             else
 334                we_patch_stag_e = we_patch_e + 1
 335             end if
 336             if (my_y /= nproc_y - 1) then
 337                sn_patch_e = sn_patch_e + 1
 338                sn_patch_stag_e = sn_patch_e
 339             else
 340                sn_patch_stag_e = sn_patch_e + 1
 341             end if
 342          else if (gridtype == 'E') then
 343             we_patch_e = we_patch_e + 1
 344             sn_patch_e = sn_patch_e + 1
 345             we_patch_stag_e = we_patch_e
 346             sn_patch_stag_e = sn_patch_e
 347          end if
 348
 349       end if
 350
 351       ! Compute multipliers for halo width; these must be 0/1
 352       if (my_x /= 0) then
 353         lh_mult = 1
 354       else
 355         lh_mult = 0
 356       end if
 357       if (my_x /= (nproc_x-1)) then
 358         rh_mult = 1
 359       else
 360         rh_mult = 0
 361       end if
 362       if (my_y /= 0) then
 363         bh_mult = 1
 364       else
 365         bh_mult = 0
 366       end if
 367       if (my_y /= (nproc_y-1)) then
 368         th_mult = 1
 369       else
 370         th_mult = 0
 371       end if
 372
 373       we_mem_s = we_patch_s - HALO_WIDTH*lh_mult
 374       we_mem_e = we_patch_e + HALO_WIDTH*rh_mult
 375       sn_mem_s = sn_patch_s - HALO_WIDTH*bh_mult
 376       sn_mem_e = sn_patch_e + HALO_WIDTH*th_mult
 377       we_mem_stag_s = we_patch_stag_s - HALO_WIDTH*lh_mult
 378       we_mem_stag_e = we_patch_stag_e + HALO_WIDTH*rh_mult
 379       sn_mem_stag_s = sn_patch_stag_s - HALO_WIDTH*bh_mult
 380       sn_mem_stag_e = sn_patch_stag_e + HALO_WIDTH*th_mult
 381
 382       ! Initialize a proj_info type for the destination grid projection. This will
 383       !   primarily be used for rotating Earth-relative winds to grid-relative winds
 384       call set_domain_projection(map_proj, stand_lon, truelat1, truelat2, &
 385                                  dom_dx, dom_dy, dom_dx, dom_dy, west_east_dim, &
 386                                  south_north_dim, real(west_east_dim)/2., &
 387                                  real(south_north_dim)/2.,cen_lat, cen_lon, &
 388                                  cen_lat, cen_lon)
 389
 390       ! Read static fields using the input module; we know that there are no more
 391       !   fields to be read when read_next_field() returns a non-zero status.
 392       istatus = 0
 393       do while (istatus == 0)
 394         call read_next_field(sp1, ep1, sp2, ep2, sp3, ep3, cname, cunits, cdesc, &
 395                              memorder, stagger, dimnames, subx, suby, &
 396                              real_array, istatus)
 397         if (istatus == 0) then
 398
 399           call mprintf(.true.,LOGFILE, 'Read in static field %s.',s1=cname)
 400           call list_insert(static_list, ckey=cname, cvalue=cname)
 401
 402           ! We will also keep copies in core of the lat/lon arrays, for use in
 403           !    interpolation of the met fields to the model grid.
 404           ! For now, we assume that the lat/lon arrays will have known field names
 405           if (index(cname, 'XLAT_M') /= 0 .and. &
 406               len_trim(cname) == len_trim('XLAT_M')) then
 407              allocate(xlat(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
 408              xlat(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
 409              call exchange_halo_r(xlat, &
 410                                   we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
 411                                   we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
 412
 413           else if (index(cname, 'XLONG_M') /= 0 .and. &
 414                    len_trim(cname) == len_trim('XLONG_M')) then
 415              allocate(xlon(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
 416              xlon(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
 417              call exchange_halo_r(xlon, &
 418                                   we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
 419                                   we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
 420
 421           else if (index(cname, 'XLAT_U') /= 0 .and. &
 422                    len_trim(cname) == len_trim('XLAT_U')) then
 423              allocate(xlat_u(we_mem_stag_s:we_mem_stag_e,sn_mem_s:sn_mem_e))
 424              xlat_u(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
 425              call exchange_halo_r(xlat_u, &
 426                                   we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
 427                                   we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
 428
 429           else if (index(cname, 'XLONG_U') /= 0 .and. &
 430                    len_trim(cname) == len_trim('XLONG_U')) then
 431              allocate(xlon_u(we_mem_stag_s:we_mem_stag_e,sn_mem_s:sn_mem_e))
 432              xlon_u(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
 433              call exchange_halo_r(xlon_u, &
 434                                   we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
 435                                   we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
 436
 437           else if (index(cname, 'XLAT_V') /= 0 .and. &
 438                    len_trim(cname) == len_trim('XLAT_V')) then
 439              allocate(xlat_v(we_mem_s:we_mem_e,sn_mem_stag_s:sn_mem_stag_e))
 440              xlat_v(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(:,:,1)
 441              call exchange_halo_r(xlat_v, &
 442                                   we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
 443                                   we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
 444
 445           else if (index(cname, 'XLONG_V') /= 0 .and. &
 446                    len_trim(cname) == len_trim('XLONG_V')) then
 447              allocate(xlon_v(we_mem_s:we_mem_e,sn_mem_stag_s:sn_mem_stag_e))
 448              xlon_v(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(:,:,1)
 449              call exchange_halo_r(xlon_v, &
 450                                   we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
 451                                   we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
 452
 453           else if (index(cname, 'LANDMASK') /= 0 .and. &
 454                    len_trim(cname) == len_trim('LANDMASK')) then
 455              allocate(landmask(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
 456              landmask(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
 457              call exchange_halo_r(landmask, &
 458                                   we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
 459                                   we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
 460
 461           end if
 462
 463           if (subx > 1) then
 464              we_mem_subgrid_s   = (we_mem_s                 + HALO_WIDTH*lh_mult - 1)*subx - HALO_WIDTH*lh_mult + 1
 465              we_mem_subgrid_e   = (we_mem_e   + (1-rh_mult) - HALO_WIDTH*rh_mult    )*subx + HALO_WIDTH*rh_mult
 466              we_patch_subgrid_s = (we_patch_s                                    - 1)*subx                      + 1
 467              we_patch_subgrid_e = (we_patch_e + (1-rh_mult)                         )*subx
 468           end if
 469           if (suby > 1) then
 470              sn_mem_subgrid_s   = (sn_mem_s                 + HALO_WIDTH*bh_mult - 1)*suby - HALO_WIDTH*bh_mult + 1
 471              sn_mem_subgrid_e   = (sn_mem_e   + (1-th_mult) - HALO_WIDTH*th_mult    )*suby + HALO_WIDTH*th_mult
 472              sn_patch_subgrid_s = (sn_patch_s                                    - 1)*suby                      + 1
 473              sn_patch_subgrid_e = (sn_patch_e + (1-th_mult)                         )*suby
 474           end if
 475
 476           ! Having read in a field, we write each level individually to the
 477           !   storage module; levels will be reassembled later on when they
 478           !   are written.
 479           do k=sp3,ep3
 480              field%header%sr_x=subx
 481              field%header%sr_y=suby
 482              field%header%version = 1
 483              field%header%date = start_date(n)
 484              field%header%time_dependent = .false.
 485              field%header%mask_field = .false.
 486              field%header%constant_field = .false.
 487              field%header%forecast_hour = 0.0
 488              field%header%fg_source = 'geogrid_model'
 489              field%header%field = cname
 490              field%header%units = cunits
 491              field%header%description = cdesc
 492              field%header%vertical_coord = dimnames(3)
 493              field%header%vertical_level = k
 494              field%header%array_order = memorder
 495              field%header%is_wind_grid_rel = .true.
 496              field%header%array_has_missing_values = .false.
 497              if (gridtype == 'C') then
 498                 if (subx > 1 .or. suby > 1) then
 499                    field%map%stagger = M
 500                    field%header%dim1(1) = we_mem_subgrid_s
 501                    field%header%dim1(2) = we_mem_subgrid_e
 502                    field%header%dim2(1) = sn_mem_subgrid_s
 503                    field%header%dim2(2) = sn_mem_subgrid_e
 504                 else if (trim(stagger) == 'M') then
 505                    field%map%stagger = M
 506                    field%header%dim1(1) = we_mem_s
 507                    field%header%dim1(2) = we_mem_e
 508                    field%header%dim2(1) = sn_mem_s
 509                    field%header%dim2(2) = sn_mem_e
 510                 else if (trim(stagger) == 'U') then
 511                    field%map%stagger = U
 512                    field%header%dim1(1) = we_mem_stag_s
 513                    field%header%dim1(2) = we_mem_stag_e
 514                    field%header%dim2(1) = sn_mem_s
 515                    field%header%dim2(2) = sn_mem_e
 516                 else if (trim(stagger) == 'V') then
 517                    field%map%stagger = V
 518                    field%header%dim1(1) = we_mem_s
 519                    field%header%dim1(2) = we_mem_e
 520                    field%header%dim2(1) = sn_mem_stag_s
 521                    field%header%dim2(2) = sn_mem_stag_e
 522                 else if (trim(stagger) == 'CORNER') then
 523                    field%map%stagger = CORNER
 524                    field%header%dim1(1) = we_mem_stag_s
 525                    field%header%dim1(2) = we_mem_stag_e
 526                    field%header%dim2(1) = sn_mem_stag_s
 527                    field%header%dim2(2) = sn_mem_stag_e
 528                 end if
 529              else if (gridtype == 'E') then
 530                 if (trim(stagger) == 'M') then
 531                    field%map%stagger = HH
 532                 else if (trim(stagger) == 'V') then
 533                    field%map%stagger = VV
 534                 end if
 535                 field%header%dim1(1) = we_mem_s
 536                 field%header%dim1(2) = we_mem_e
 537                 field%header%dim2(1) = sn_mem_s
 538                 field%header%dim2(2) = sn_mem_e
 539              end if
 540
 541              allocate(field%valid_mask)
 542
 543              if (subx > 1 .or. suby > 1) then
 544                 allocate(field%r_arr(we_mem_subgrid_s:we_mem_subgrid_e,&
 545                                      sn_mem_subgrid_s:sn_mem_subgrid_e))
 546                 field%r_arr(we_patch_subgrid_s:we_patch_subgrid_e,sn_patch_subgrid_s:sn_patch_subgrid_e) = &
 547                            real_array(sp1:ep1,sp2:ep2,k)
 548                 call exchange_halo_r(field%r_arr, &
 549                            we_mem_subgrid_s, we_mem_subgrid_e, sn_mem_subgrid_s, sn_mem_subgrid_e, 1, 1, &
 550                            we_patch_subgrid_s, we_patch_subgrid_e, sn_patch_subgrid_s, sn_patch_subgrid_e, 1, 1)
 551                 call bitarray_create(field%valid_mask, &
 552                                      (we_mem_subgrid_e-we_mem_subgrid_s)+1, &
 553                                      (sn_mem_subgrid_e-sn_mem_subgrid_s)+1)
 554                 do j=1,(sn_mem_subgrid_e-sn_mem_subgrid_s)+1
 555                    do i=1,(we_mem_subgrid_e-we_mem_subgrid_s)+1
 556                       call bitarray_set(field%valid_mask, i, j)
 557                    end do
 558                 end do
 559
 560              else if (field%map%stagger == M  .or. &
 561                  field%map%stagger == HH .or. &
 562                  field%map%stagger == VV) then
 563                 allocate(field%r_arr(we_mem_s:we_mem_e,&
 564                                      sn_mem_s:sn_mem_e))
 565                 field%r_arr(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(sp1:ep1,sp2:ep2,k)
 566                 call exchange_halo_r(field%r_arr, &
 567                            we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
 568                            we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
 569                 call bitarray_create(field%valid_mask, &
 570                                      (we_mem_e-we_mem_s)+1, &
 571                                      (sn_mem_e-sn_mem_s)+1)
 572                 do j=1,(sn_mem_e-sn_mem_s)+1
 573                    do i=1,(we_mem_e-we_mem_s)+1
 574                       call bitarray_set(field%valid_mask, i, j)
 575                    end do
 576                 end do
 577              else if (field%map%stagger == U) then
 578                 allocate(field%r_arr(we_mem_stag_s:we_mem_stag_e,&
 579                                      sn_mem_s:sn_mem_e))
 580                 field%r_arr(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(sp1:ep1,sp2:ep2,k)
 581                 call exchange_halo_r(field%r_arr, &
 582                            we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
 583                            we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
 584                 call bitarray_create(field%valid_mask, &
 585                                      (we_mem_stag_e-we_mem_stag_s)+1, &
 586                                      (sn_mem_e-sn_mem_s)+1)
 587                 do j=1,(sn_mem_e-sn_mem_s)+1
 588                    do i=1,(we_mem_stag_e-we_mem_stag_s)+1
 589                       call bitarray_set(field%valid_mask, i, j)
 590                    end do
 591                 end do
 592              else if (field%map%stagger == V) then
 593                 allocate(field%r_arr(we_mem_s:we_mem_e,&
 594                                      sn_mem_stag_s:sn_mem_stag_e))
 595                 field%r_arr(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(sp1:ep1,sp2:ep2,k)
 596                 call exchange_halo_r(field%r_arr, &
 597                            we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
 598                            we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
 599                 call bitarray_create(field%valid_mask, &
 600                                      (we_mem_e-we_mem_s)+1, &
 601                                      (sn_mem_stag_e-sn_mem_stag_s)+1)
 602                 do j=1,(sn_mem_stag_e-sn_mem_stag_s)+1
 603                    do i=1,(we_mem_e-we_mem_s)+1
 604                       call bitarray_set(field%valid_mask, i, j)
 605                    end do
 606                 end do
 607              else if (field%map%stagger == CORNER) then
 608                 allocate(field%r_arr(we_mem_stag_s:we_mem_stag_e,&
 609                                      sn_mem_stag_s:sn_mem_stag_e))
 610                 field%r_arr(we_patch_stag_s:we_patch_stag_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(sp1:ep1,sp2:ep2,k)
 611                 call exchange_halo_r(field%r_arr, &
 612                            we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
 613                            we_patch_stag_s, we_patch_stag_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
 614                 call bitarray_create(field%valid_mask, &
 615                                      (we_mem_stag_e-we_mem_stag_s)+1, &
 616                                      (sn_mem_stag_e-sn_mem_stag_s)+1)
 617                 do j=1,(sn_mem_stag_e-sn_mem_stag_s)+1
 618                    do i=1,(we_mem_stag_e-we_mem_stag_s)+1
 619                       call bitarray_set(field%valid_mask, i, j)
 620                    end do
 621                 end do
 622              end if
 623
 624              nullify(field%modified_mask)
 625
 626              call storage_put_field(field)
 627
 628           end do
 629
 630         end if   ! if (istatus == 0)
 631       end do  ! do while (istatus == 0)
 632
 633       static_list_array => list_get_keys(static_list)
 634       call list_init(flag_list)
 635       do i=1,size(static_list_array)
 636          istatus = 0
 637          ctemp = 'FLAG_'//trim(static_list_array(i)%ckey)
 638          call ext_get_dom_ti_integer_scalar(ctemp, istatus, suppress_errors=.true.)
 639          if (istatus == 1) call list_insert(flag_list, ckey=ctemp, cvalue=ctemp)
 640       end do
 641       deallocate(static_list_array)
 642       call list_destroy(static_list)
 643
 644       ! Done reading all static fields for this domain
 645       call input_close()
 646
 647       static_list_array => list_get_keys(flag_list)
 648       allocate(geogrid_flags(size(static_list_array)))
 649       do i=1,size(static_list_array)
 650          geogrid_flags(i) = static_list_array(i)%ckey
 651       end do
 652       deallocate(static_list_array)
 653
 654       call list_destroy(flag_list)
 655
 656    end subroutine get_static_fields
 657
 658
 659    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 660    ! Name: process_single_met_time
 661    !
 662    ! Purpose:
 663    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 664    subroutine process_single_met_time(do_const_processing, &
 665                              temp_date, n, extra_row, extra_col, xlat, xlon, &
 666                              xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
 667                              title, dyn_opt, &
 668                              west_east_dim, south_north_dim, &
 669                              we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
 670                              we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
 671                              sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
 672                              we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
 673                              sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
 674                              got_this_field, &
 675                              map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
 676                              grid_id, parent_id, i_parent_start, &
 677                              j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
 678                              cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
 679                              pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
 680                              corner_lats, corner_lons, output_flags, geogrid_flags, process_bdy_width)
 681
 682       use bitarray_module
 683       use gridinfo_module
 684       use interp_module
 685       use interp_option_module
 686       use llxy_module
 687       use misc_definitions_module
 688       use module_debug
 689       use output_module
 690       use parallel_module
 691       use read_met_module
 692       use rotate_winds_module
 693       use storage_module
 694
 695       implicit none
 696
 697       ! Arguments
 698       logical, intent(in) :: do_const_processing
 699       integer, intent(in) :: n, dyn_opt, west_east_dim, south_north_dim, map_proj, &
 700                  we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
 701                  we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
 702                  sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
 703                  we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
 704                  sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
 705                  is_water, is_lake, is_ice, is_urban, i_soilwater, &
 706                  grid_id, parent_id, i_parent_start, j_parent_start, &
 707                  i_parent_end, j_parent_end, parent_grid_ratio, sub_x, sub_y, num_land_cat, &
 708                  process_bdy_width
 709 ! BUG: Should we be passing these around as pointers, or just declare them as arrays?
 710       real, pointer, dimension(:,:) :: landmask
 711       real, intent(in) :: dom_dx, dom_dy, cen_lat, moad_cen_lat, cen_lon, stand_lon, &
 712                  truelat1, truelat2, pole_lat, pole_lon
 713       real, dimension(16), intent(in) :: corner_lats, corner_lons
 714       real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
 715       logical, intent(in) :: extra_row, extra_col
 716       logical, dimension(:), intent(inout) :: got_this_field
 717       character (len=19), intent(in) :: temp_date
 718       character (len=128), intent(in) :: mminlu
 719       character (len=128), dimension(:), intent(inout) :: output_flags
 720       character (len=128), dimension(:), pointer :: geogrid_flags
 721
 722 ! BUG: Move this constant to misc_definitions_module?
 723 integer, parameter :: BDR_WIDTH = 3
 724
 725       ! Local variables
 726       integer :: istatus, fg_idx, idx, idxt, i, j, bottom_top_dim, &
 727                  sm1, em1, sm2, em2, sm3, em3, &
 728                  sp1, ep1, sp2, ep2, sp3, ep3, &
 729                  sd1, ed1, sd2, ed2, sd3, ed3, &
 730                  u_idx
 731       integer :: nmet_flags
 732       integer :: num_metgrid_soil_levs
 733       integer, pointer, dimension(:) :: soil_levels
 734       real :: rx, ry
 735       integer, pointer, dimension(:) :: u_levels, v_levels
 736       real, pointer, dimension(:,:,:) :: real_array
 737       character (len=19) :: output_date
 738       character (len=128) :: cname, title
 739       character (len=MAX_FILENAME_LEN) :: input_name
 740       character (len=128), allocatable, dimension(:) :: met_flags
 741       type (fg_input) :: field, u_field, v_field
 742       type (met_data) :: fg_data
 743
 744       ! CWH Initialize local pointer variables
 745       nullify(soil_levels)
 746       nullify(u_levels)
 747       nullify(v_levels)
 748       nullify(real_array)
 749
 750
 751       ! For this time, we need to process all first-guess filename roots. When we
 752       !   hit a root containing a '*', we assume we have hit the end of the list
 753       fg_idx = 1
 754       if (do_const_processing) then
 755          input_name = constants_name(fg_idx)
 756       else
 757          input_name = fg_name(fg_idx)
 758       end if
 759       do while (input_name /= '*')
 760
 761          call mprintf(.true.,STDOUT, '    %s', s1=input_name)
 762          call mprintf(.true.,LOGFILE, 'Getting input fields from %s', s1=input_name)
 763
 764          if (index(input_name, 'mpas:') == 1) then
 765             call process_mpas_fields(input_name, do_const_processing, temp_date, fg_data, got_this_field, &
 766                                      landmask, process_bdy_width, &
 767                                      u_field, v_field, &
 768                                      dom_dx, dom_dy, &
 769                                      xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v, &
 770                                      output_flags, west_east_dim, south_north_dim, &
 771                                      we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
 772                                      we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
 773                                      sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
 774                                      we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
 775                                      sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e )
 776          else
 777             call process_intermediate_fields(input_name, do_const_processing, temp_date, fg_data, got_this_field, &
 778                                              landmask, process_bdy_width, &
 779                                              u_field, v_field, &
 780                                              dom_dx, dom_dy, &
 781                                              xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v, &
 782                                              output_flags, west_east_dim, south_north_dim, &
 783                                              we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
 784                                              we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
 785                                              sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
 786                                              we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
 787                                              sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e )
 788          end if
 789
 790          fg_idx = fg_idx + 1
 791          if (do_const_processing) then
 792             input_name = constants_name(fg_idx)
 793          else
 794             input_name = fg_name(fg_idx)
 795          end if
 796       end do ! while (input_name /= '*')
 797
 798       !
 799       ! Rotate winds from earth-relative to grid-relative
 800       !
 801
 802       call storage_get_levels(u_field, u_levels)
 803       call storage_get_levels(v_field, v_levels)
 804
 805       if (associated(u_levels) .and. associated(v_levels)) then
 806          u_idx = 1
 807          do u_idx = 1, size(u_levels)
 808             u_field%header%vertical_level = u_levels(u_idx)
 809             call storage_get_field(u_field, istatus)
 810             v_field%header%vertical_level = v_levels(u_idx)
 811             call storage_get_field(v_field, istatus)
 812
 813             if (gridtype == 'C') then
 814                call met_to_map(u_field%r_arr, u_field%valid_mask, &
 815                                v_field%r_arr, v_field%valid_mask, &
 816                                we_mem_stag_s, sn_mem_s, &
 817                                we_mem_stag_e, sn_mem_e, &
 818                                we_mem_s, sn_mem_stag_s, &
 819                                we_mem_e, sn_mem_stag_e, &
 820                                xlon_u, xlon_v, xlat_u, xlat_v)
 821             else if (gridtype == 'E') then
 822                call met_to_map_nmm(u_field%r_arr, u_field%valid_mask, &
 823                                    v_field%r_arr, v_field%valid_mask, &
 824                                    we_mem_s, sn_mem_s, &
 825                                    we_mem_e, sn_mem_e, &
 826                                    xlat_v, xlon_v)
 827             end if
 828
 829          end do
 830
 831          deallocate(u_levels)
 832          deallocate(v_levels)
 833
 834       end if
 835
 836       if (do_const_processing) return
 837
 838       !
 839       ! Now that we have all degribbed fields, we build a 3-d pressure field, and fill in any
 840       !   missing levels in the other 3-d fields
 841       !
 842       call mprintf(.true.,LOGFILE,'Filling missing levels.')
 843       call fill_missing_levels(output_flags)
 844
 845       call mprintf(.true.,LOGFILE,'Creating derived fields.')
 846       call create_derived_fields(gridtype, fg_data%hdate, fg_data%xfcst, &
 847                                  we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
 848                                  we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e, &
 849                                  got_this_field, output_flags)
 850
 851       !
 852       ! Derive some MPAS fields from others, e.g., TT from theta and pressure
 853       !
 854       call derive_mpas_fields()
 855
 856       !
 857       ! Check that every mandatory field was found in input data
 858       !
 859       do i=1,num_entries
 860          if (is_mandatory(i) .and. .not. got_this_field(i)) then
 861             call mprintf(.true.,ERROR,'The mandatory field %s was not found in any input data.',s1=fieldname(i))
 862          end if
 863       end do
 864
 865       !
 866       ! Before we begin to write fields, if debug_level is set high enough, we
 867       !    write a table of which fields are available at which levels to the
 868       !    metgrid.log file, and then we check to see if any fields are not
 869       !    completely covered with data.
 870       !
 871       call storage_print_fields()
 872       call find_missing_values()
 873
 874       !
 875       ! All of the processing is now done for this time period for this domain;
 876       !   now we simply output every field from the storage module.
 877       !
 878
 879       title = 'OUTPUT FROM METGRID V4.4'
 880
 881       ! Initialize the output module for this domain and time
 882       call mprintf(.true.,LOGFILE,'Initializing output module.')
 883       output_date = temp_date
 884       if ( .not. nocolons ) then
 885          if (len_trim(temp_date) == 13) then
 886             output_date(14:19) = ':00:00'
 887          else if (len_trim(temp_date) == 16) then
 888             output_date(17:19) = ':00'
 889          end if
 890       else
 891          if (len_trim(temp_date) == 13) then
 892             output_date(14:19) = '_00_00'
 893          else if (len_trim(temp_date) == 16) then
 894             output_date(17:19) = '_00'
 895          end if
 896       endif
 897
 898       call output_init(n, title, output_date, gridtype, dyn_opt, &
 899                        corner_lats, corner_lons, &
 900                        we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
 901                        we_patch_s,  we_patch_e,  sn_patch_s,  sn_patch_e, &
 902                        we_mem_s,    we_mem_e,    sn_mem_s,    sn_mem_e, &
 903                        extra_col, extra_row)
 904
 905       call get_bottom_top_dim(bottom_top_dim)
 906
 907       ! Add in a flag to tell real that we have seven new msf fields
 908       nmet_flags = num_entries + 1
 909       if (associated(geogrid_flags)) nmet_flags = nmet_flags + size(geogrid_flags)
 910       allocate(met_flags(nmet_flags))
 911       do i=1,num_entries
 912          met_flags(i) = output_flags(i)
 913       end do
 914       if (gridtype == 'C') then
 915          met_flags(num_entries+1) = 'FLAG_MF_XY'
 916       else
 917          met_flags(num_entries+1) = ' '
 918       end if
 919       if (associated(geogrid_flags)) then
 920          do i=1,size(geogrid_flags)
 921             met_flags(num_entries+1+i) = geogrid_flags(i)
 922          end do
 923       end if
 924
 925       ! Find out how many soil levels or layers we have; this assumes a field named ST
 926       field % header % field = 'ST'
 927       nullify(soil_levels)
 928       call storage_get_levels(field, soil_levels)
 929
 930       if (.not. associated(soil_levels)) then
 931          field % header % field = 'SOILT'
 932          nullify(soil_levels)
 933          call storage_get_levels(field, soil_levels)
 934       end if
 935
 936       if (.not. associated(soil_levels)) then
 937          field % header % field = 'STC_WPS'
 938          nullify(soil_levels)
 939          call storage_get_levels(field, soil_levels)
 940       end if
 941
 942       if (associated(soil_levels)) then
 943          num_metgrid_soil_levs = size(soil_levels)
 944          deallocate(soil_levels)
 945       else
 946          num_metgrid_soil_levs = 0
 947       end if
 948
 949       ! First write out global attributes
 950       call mprintf(.true.,LOGFILE,'Writing global attributes to output.')
 951       call write_global_attrs(title, output_date, gridtype, dyn_opt, west_east_dim,          &
 952                               south_north_dim, bottom_top_dim,                               &
 953                               we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e,      &
 954                               sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e,      &
 955                               map_proj, mminlu, num_land_cat,                                &
 956                               is_water, is_lake, is_ice, is_urban, i_soilwater,              &
 957                               grid_id, parent_id, i_parent_start,                            &
 958                               j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy,    &
 959                               cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
 960                               pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y,           &
 961                               corner_lats, corner_lons, num_metgrid_soil_levs,               &
 962                               met_flags, nmet_flags, process_bdy_width)
 963
 964       deallocate(met_flags)
 965
 966       call reset_next_field()
 967
 968       istatus = 0
 969
 970       ! Now loop over all output fields, writing each to the output module
 971       do while (istatus == 0)
 972          call get_next_output_field(cname, real_array, &
 973                                     sm1, em1, sm2, em2, sm3, em3, istatus)
 974          if (istatus == 0) then
 975
 976             call mprintf(.true.,LOGFILE,'Writing field %s to output.',s1=cname)
 977             call write_field(sm1, em1, sm2, em2, sm3, em3, &
 978                              cname, output_date, real_array)
 979             deallocate(real_array)
 980
 981          end if
 982
 983       end do
 984
 985       call mprintf(.true.,LOGFILE,'Closing output file.')
 986       call output_close()
 987
 988       ! Free up memory used by met fields for this valid time
 989       call storage_delete_all_td()
 990
 991    end subroutine process_single_met_time
 992
 993
 994    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 995    ! Name: process_intermediate_fields
 996    !
 997    ! Purpose:
 998    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 999    subroutine process_intermediate_fields(input_name, do_const_processing, temp_date, fg_data, got_this_field, &
1000                                           landmask, process_bdy_width, &
1001                                           u_field, v_field, &
1002                                           dom_dx, dom_dy, &
1003                                           xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v, &
1004                                           output_flags, west_east_dim, south_north_dim, &
1005                                           we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1006                                           we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
1007                                           sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
1008                                           we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
1009                                           sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e )
1010
1011       use bitarray_module
1012       use gridinfo_module
1013       use interp_module
1014       use interp_option_module
1015       use llxy_module
1016       use misc_definitions_module
1017       use module_debug
1018       use output_module
1019       use parallel_module
1020       use read_met_module
1021       use rotate_winds_module
1022       use storage_module
1023
1024       implicit none
1025
1026 ! BUG: Move this constant to misc_definitions_module?
1027 integer, parameter :: BDR_WIDTH = 3
1028
1029       character (len=*), intent(inout) :: input_name
1030       logical, intent(in) :: do_const_processing
1031       character (len=*), intent(in) :: temp_date
1032       type (met_data), intent(inout) :: fg_data
1033       logical, dimension(:), intent(inout) :: got_this_field
1034       real, pointer, dimension(:,:) :: landmask
1035       integer, intent(in) :: process_bdy_width
1036       type (fg_input), intent(inout) :: u_field, v_field
1037       character (len=128), dimension(:), intent(inout) :: output_flags
1038       real, intent(in) :: dom_dx, dom_dy
1039       real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
1040       integer, intent(in) :: west_east_dim, south_north_dim, &
1041                  we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1042                  we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
1043                  sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
1044                  we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
1045                  sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e
1046
1047       integer :: istatus
1048       integer :: idx, idxt, u_idx
1049       integer :: iqstatus
1050       real :: threshold
1051       real, pointer, dimension(:,:) :: halo_slab => null()
1052       integer, pointer, dimension(:) :: u_levels => null()
1053       integer, pointer, dimension(:) :: v_levels => null()
1054       integer :: bdr_wdth
1055       logical :: do_gcell_interp
1056       type (fg_input) :: field
1057
1058
1059       ! Do a first pass through this fg source to get all mask fields used
1060       !   during interpolation
1061       call get_interp_masks(trim(input_name), do_const_processing, temp_date)
1062
1063       istatus = 0
1064
1065       ! Initialize the module for reading in the met fields
1066       call read_met_init(trim(input_name), do_const_processing, temp_date, istatus)
1067
1068       if (istatus == 0) then
1069
1070          ! Process all fields and levels from the current file; read_next_met_field()
1071          !   will return a non-zero status when there are no more fields to be read.
1072          do while (istatus == 0)
1073
1074             call read_next_met_field(fg_data, istatus)
1075
1076             if (istatus == 0) then
1077
1078                ! Find index into fieldname, interp_method, masked, and fill_missing
1079                !   of the current field
1080                idxt = num_entries + 1
1081                do idx=1,num_entries
1082                   if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
1083                       (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
1084
1085                      got_this_field(idx) = .true.
1086
1087                      if (index(input_name,trim(from_input(idx))) /= 0 .or. &
1088                         (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1089                         idxt = idx
1090                      end if
1091
1092                   end if
1093                end do
1094                idx = idxt
1095                if (idx > num_entries) idx = num_entries ! The last entry is a default
1096
1097                ! Do we need to rename this field?
1098                if (output_name(idx) /= ' ') then
1099                   fg_data%field = output_name(idx)(1:9)
1100
1101                   idxt = num_entries + 1
1102                   do idx=1,num_entries
1103                      if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
1104                          (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
1105
1106                         got_this_field(idx) = .true.
1107
1108                         if (index(input_name,trim(from_input(idx))) /= 0 .or. &
1109                            (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1110                            idxt = idx
1111                         end if
1112
1113                      end if
1114                   end do
1115                   idx = idxt
1116                   if (idx > num_entries) idx = num_entries ! The last entry is a default
1117                end if
1118
1119                ! Do a simple check to see whether this is a global dataset
1120                ! Note that we do not currently support regional Gaussian grids
1121                if ((fg_data%iproj == PROJ_LATLON .and. abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) &
1122                     .or. (fg_data%iproj == PROJ_GAUSS)) then
1123                   bdr_wdth = BDR_WIDTH
1124                   allocate(halo_slab(1-BDR_WIDTH:fg_data%nx+BDR_WIDTH,1:fg_data%ny))
1125
1126                   halo_slab(1:fg_data%nx,                      1:fg_data%ny) = &
1127                             fg_data%slab(1:fg_data%nx,              1:fg_data%ny)
1128
1129                   halo_slab(1-BDR_WIDTH:0,                     1:fg_data%ny) = &
1130                             fg_data%slab(fg_data%nx-BDR_WIDTH+1:fg_data%nx, 1:fg_data%ny)
1131
1132                   halo_slab(fg_data%nx+1:fg_data%nx+BDR_WIDTH, 1:fg_data%ny) = &
1133                             fg_data%slab(1:BDR_WIDTH,       1:fg_data%ny)
1134
1135                   deallocate(fg_data%slab)
1136                else
1137                   bdr_wdth = 0
1138                   halo_slab => fg_data%slab
1139                   nullify(fg_data%slab)
1140                end if
1141
1142                call mprintf(.true.,LOGFILE,'Processing %s at level %f.',s1=fg_data%field,f1=fg_data%xlvl)
1143
1144                call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, &
1145                                            fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, &
1146                                            fg_data%deltalon, fg_data%starti, fg_data%startj, &
1147                                            fg_data%startlat, fg_data%startlon,  &
1148                                            fg_data%pole_lat, fg_data%pole_lon,        &
1149                                            fg_data%centerlat, fg_data%centerlon,      &
1150                                            real(fg_data%nx+1)/2., real(fg_data%ny+1)/2.,  &
1151                                            earth_radius=fg_data%earth_radius*1000.)
1152
1153                ! Initialize fg_input structure to store the field
1154                field%header%version = 1
1155                field%header%date = fg_data%hdate//'        '
1156                if (do_const_processing) then
1157                   field%header%time_dependent = .false.
1158                   field%header%constant_field = .true.
1159                else
1160                   field%header%time_dependent = .true.
1161                   field%header%constant_field = .false.
1162                end if
1163                field%header%forecast_hour = fg_data%xfcst
1164                field%header%fg_source = 'FG'
1165                field%header%field = ' '
1166                field%header%field(1:9) = fg_data%field
1167                field%header%units = ' '
1168                field%header%units(1:25) = fg_data%units
1169                field%header%description = ' '
1170                field%header%description(1:46) = fg_data%desc
1171                call get_z_dim_name(fg_data%field,field%header%vertical_coord)
1172                field%header%vertical_level = nint(fg_data%xlvl)
1173                field%header%sr_x = 1
1174                field%header%sr_y = 1
1175                field%header%array_order = 'XY '
1176                field%header%is_wind_grid_rel = fg_data%is_wind_grid_rel
1177                field%header%array_has_missing_values = .false.
1178                nullify(field%r_arr)
1179                nullify(field%valid_mask)
1180                nullify(field%modified_mask)
1181
1182                if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
1183                   output_flags(idx) = flag_in_output(idx)
1184                end if
1185
1186                ! If we should not output this field, just list it as a mask field
1187                if (output_this_field(idx)) then
1188                   field%header%mask_field = .false.
1189                else
1190                   field%header%mask_field = .true.
1191                end if
1192
1193                !
1194                ! Before actually doing any interpolation to the model grid, we must check
1195                !    whether we will be using the average_gcell interpolator that averages all
1196                !    source points in each model grid cell
1197                !
1198                do_gcell_interp = .false.
1199                if (index(interp_method(idx),'average_gcell') /= 0) then
1200
1201                   call get_gcell_threshold(interp_method(idx), threshold, istatus)
1202                   if (istatus == 0) then
1203                      if (fg_data%dx == 0. .and. fg_data%dy == 0. .and. &
1204                          fg_data%deltalat /= 0. .and. fg_data%deltalon /= 0.) then
1205                         fg_data%dx = abs(fg_data%deltalon)
1206                         fg_data%dy = abs(fg_data%deltalat)
1207                      else
1208 ! BUG: Need to more correctly handle dx/dy in meters.
1209                         fg_data%dx = fg_data%dx / 111000.  ! Convert meters to approximate degrees
1210                         fg_data%dy = fg_data%dy / 111000.
1211                      end if
1212                      if (gridtype == 'C') then
1213                         if (threshold*max(fg_data%dx,fg_data%dy)*111. <= max(dom_dx,dom_dy)/1000.) &
1214                            do_gcell_interp = .true.
1215                      else if (gridtype == 'E') then
1216                         if (threshold*max(fg_data%dx,fg_data%dy) <= max(dom_dx,dom_dy)) &
1217                            do_gcell_interp = .true.
1218                      end if
1219                   end if
1220                end if
1221
1222                ! Interpolate to U staggering
1223                if (output_stagger(idx) == U) then
1224
1225                   call storage_query_field(field, iqstatus)
1226                   if (iqstatus == 0) then
1227                      call storage_get_field(field, iqstatus)
1228                      call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
1229                                   ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
1230                      if (associated(field%modified_mask)) then
1231                         call bitarray_destroy(field%modified_mask)
1232                         nullify(field%modified_mask)
1233                      end if
1234                   else
1235                      allocate(field%valid_mask)
1236                      call bitarray_create(field%valid_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
1237                   end if
1238
1239                   ! Save a copy of the fg_input structure for the U field so that we can find it later
1240                   if (is_u_field(idx)) call dup(field, u_field)
1241
1242                   allocate(field%modified_mask)
1243                   call bitarray_create(field%modified_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
1244
1245                   if (do_const_processing .or. field%header%time_dependent) then
1246                      call interp_met_field(input_name, fg_data%field, U, M, &
1247                                   field, xlat_u, xlon_u, we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, &
1248                                   we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1249                                   halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
1250                                   field%modified_mask, process_bdy_width)
1251                   else
1252                      call mprintf(.true.,INFORM,' - already processed this field from constant file.')
1253                   end if
1254
1255                ! Interpolate to V staggering
1256                else if (output_stagger(idx) == V) then
1257
1258                   call storage_query_field(field, iqstatus)
1259                   if (iqstatus == 0) then
1260                      call storage_get_field(field, iqstatus)
1261                      call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
1262                                   ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
1263                      if (associated(field%modified_mask)) then
1264                         call bitarray_destroy(field%modified_mask)
1265                         nullify(field%modified_mask)
1266                      end if
1267                   else
1268                      allocate(field%valid_mask)
1269                      call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
1270                   end if
1271
1272                   ! Save a copy of the fg_input structure for the V field so that we can find it later
1273                   if (is_v_field(idx)) call dup(field, v_field)
1274
1275                   allocate(field%modified_mask)
1276                   call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
1277
1278                   if (do_const_processing .or. field%header%time_dependent) then
1279                      call interp_met_field(input_name, fg_data%field, V, M, &
1280                                   field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
1281                                   we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1282                                   halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
1283                                   field%modified_mask, process_bdy_width)
1284                   else
1285                      call mprintf(.true.,INFORM,' - already processed this field from constant file.')
1286                   end if
1287
1288                ! Interpolate to VV staggering
1289                else if (output_stagger(idx) == VV) then
1290
1291                   call storage_query_field(field, iqstatus)
1292                   if (iqstatus == 0) then
1293                      call storage_get_field(field, iqstatus)
1294                      call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
1295                                   ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
1296                      if (associated(field%modified_mask)) then
1297                         call bitarray_destroy(field%modified_mask)
1298                         nullify(field%modified_mask)
1299                      end if
1300                   else
1301                      allocate(field%valid_mask)
1302                      call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
1303                   end if
1304
1305                   ! Save a copy of the fg_input structure for the U field so that we can find it later
1306                   if (is_u_field(idx)) call dup(field, u_field)
1307
1308                   ! Save a copy of the fg_input structure for the V field so that we can find it later
1309                   if (is_v_field(idx)) call dup(field, v_field)
1310
1311                   allocate(field%modified_mask)
1312                   call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
1313
1314                   if (do_const_processing .or. field%header%time_dependent) then
1315                      call interp_met_field(input_name, fg_data%field, VV, M, &
1316                                   field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1317                                   we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1318                                   halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
1319                                   field%modified_mask, process_bdy_width)
1320                   else
1321                      call mprintf(.true.,INFORM,' - already processed this field from constant file.')
1322                   end if
1323
1324                ! All other fields interpolated to M staggering for C grid, H staggering for E grid
1325                else
1326
1327                   call storage_query_field(field, iqstatus)
1328                   if (iqstatus == 0) then
1329                      call storage_get_field(field, iqstatus)
1330                      call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
1331                                   ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
1332                      if (associated(field%modified_mask)) then
1333                         call bitarray_destroy(field%modified_mask)
1334                         nullify(field%modified_mask)
1335                      end if
1336                   else
1337                      allocate(field%valid_mask)
1338                      call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
1339                   end if
1340
1341                   allocate(field%modified_mask)
1342                   call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
1343
1344                   if (do_const_processing .or. field%header%time_dependent) then
1345                      if (gridtype == 'C') then
1346                         call interp_met_field(input_name, fg_data%field, M, M, &
1347                                      field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1348                                      we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1349                                      halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
1350                                      field%modified_mask, process_bdy_width, landmask)
1351
1352                      else if (gridtype == 'E') then
1353                         call interp_met_field(input_name, fg_data%field, HH, M, &
1354                                      field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1355                                      we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1356                                      halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
1357                                      field%modified_mask, process_bdy_width, landmask)
1358                      end if
1359                   else
1360                      call mprintf(.true.,INFORM,' - already processed this field from constant file.')
1361                   end if
1362
1363                end if
1364
1365                call bitarray_merge(field%valid_mask, field%modified_mask)
1366
1367                deallocate(halo_slab)
1368
1369                ! Store the interpolated field
1370                call storage_put_field(field)
1371
1372                call pop_source_projection()
1373
1374             end if
1375          end do
1376
1377          call read_met_close()
1378
1379          call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, &
1380                                      fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, &
1381                                      fg_data%deltalon, fg_data%starti, fg_data%startj, &
1382                                      fg_data%startlat, fg_data%startlon,  &
1383                                      fg_data%pole_lat, fg_data%pole_lon,        &
1384                                      fg_data%centerlat, fg_data%centerlon,      &
1385                                      real(fg_data%nx+1)/2., real(fg_data%ny+1)/2.,  &
1386                                      earth_radius=fg_data%earth_radius*1000.)
1387
1388          !
1389          ! If necessary, rotate winds to earth-relative for this fg source
1390          !
1391
1392          call storage_get_levels(u_field, u_levels)
1393          call storage_get_levels(v_field, v_levels)
1394
1395          if (associated(u_levels) .and. associated(v_levels)) then
1396             u_idx = 1
1397             do u_idx = 1, size(u_levels)
1398                u_field%header%vertical_level = u_levels(u_idx)
1399                call storage_get_field(u_field, istatus)
1400                v_field%header%vertical_level = v_levels(u_idx)
1401                call storage_get_field(v_field, istatus)
1402
1403                if (associated(u_field%modified_mask) .and. &
1404                    associated(v_field%modified_mask)) then
1405
1406                   if (u_field%header%is_wind_grid_rel) then
1407                      if (gridtype == 'C') then
1408                         call map_to_met(u_field%r_arr, u_field%modified_mask, &
1409                                         v_field%r_arr, v_field%modified_mask, &
1410                                         we_mem_stag_s, sn_mem_s, &
1411                                         we_mem_stag_e, sn_mem_e, &
1412                                         we_mem_s, sn_mem_stag_s, &
1413                                         we_mem_e, sn_mem_stag_e, &
1414                                         xlon_u, xlon_v, xlat_u, xlat_v)
1415                      else if (gridtype == 'E') then
1416                         call map_to_met_nmm(u_field%r_arr, u_field%modified_mask, &
1417                                             v_field%r_arr, v_field%modified_mask, &
1418                                             we_mem_s, sn_mem_s, &
1419                                             we_mem_e, sn_mem_e, &
1420                                             xlat_v, xlon_v)
1421                      end if
1422                   end if
1423
1424                   call bitarray_destroy(u_field%modified_mask)
1425                   call bitarray_destroy(v_field%modified_mask)
1426                   nullify(u_field%modified_mask)
1427                   nullify(v_field%modified_mask)
1428                   call storage_put_field(u_field)
1429                   call storage_put_field(v_field)
1430                end if
1431
1432             end do
1433
1434             deallocate(u_levels)
1435             deallocate(v_levels)
1436
1437          end if
1438
1439          call pop_source_projection()
1440
1441       else
1442          if (do_const_processing) then
1443             call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=input_name)
1444          else
1445             call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=trim(input_name)//':'//trim(temp_date))
1446          end if
1447       end if
1448
1449    end subroutine process_intermediate_fields
1450
1451
1452    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1453    ! Name: process_mpas_fields
1454    !
1455    ! Purpose:
1456    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1457    subroutine process_mpas_fields(input_name, do_const_processing, temp_date, fg_data, got_this_field, &
1458                                           landmask, process_bdy_width, &
1459                                           u_field, v_field, &
1460                                           dom_dx, dom_dy, &
1461                                           xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v, &
1462                                           output_flags, west_east_dim, south_north_dim, &
1463                                           we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1464                                           we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
1465                                           sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
1466                                           we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
1467                                           sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e )
1468
1469       use bitarray_module
1470       use gridinfo_module
1471       use interp_module
1472       use interp_option_module
1473       use read_met_module
1474       use llxy_module
1475       use misc_definitions_module
1476       use module_debug
1477       use output_module
1478       use parallel_module
1479       use rotate_winds_module
1480       use storage_module
1481       use scan_input
1482       use mpas_mesh
1483
1484       implicit none
1485
1486 ! BUG: Move this constant to misc_definitions_module?
1487 integer, parameter :: BDR_WIDTH = 3
1488
1489       character (len=*), intent(inout) :: input_name
1490       logical, intent(in) :: do_const_processing
1491       character (len=*), intent(in) :: temp_date
1492       type (met_data), intent(inout) :: fg_data
1493       logical, dimension(:), intent(inout) :: got_this_field
1494       real, pointer, dimension(:,:) :: landmask
1495       integer, intent(in) :: process_bdy_width
1496       type (fg_input), intent(inout) :: u_field, v_field
1497       character (len=128), dimension(:), intent(inout) :: output_flags
1498       real, intent(in) :: dom_dx, dom_dy
1499       real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
1500       integer, intent(in) :: west_east_dim, south_north_dim, &
1501                  we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1502                  we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
1503                  sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
1504                  we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
1505                  sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e
1506
1507       real, parameter :: deg2rad = asin(1.0) / 90.0
1508
1509       integer :: i, j, k
1510       integer :: idx
1511       integer :: istat
1512       integer :: strlen
1513       character (len=MAX_FILENAME_LEN) :: mpas_filename
1514       integer :: nRecords
1515       type (input_handle_type) :: mpas_handle
1516       type (input_field_type) :: mpas_field
1517       type (target_field_type) :: wrf_field
1518       type (fg_input) :: field_to_store
1519
1520       strlen = len_trim(input_name)
1521       if (do_const_processing) then
1522          write(mpas_filename,'(a)') input_name(6:strlen)
1523       else
1524          write(mpas_filename,'(a)') input_name(6:strlen)//'.'//trim(temp_date)//'.nc'
1525       end if
1526       call mprintf(.true.,LOGFILE,'Processing MPAS fields from file %s',s1=mpas_filename)
1527
1528       !
1529       ! If we do not already have mesh information, get that now...
1530       !
1531       if (.not. mpas_source_mesh % valid) then
1532          if (mpas_mesh_setup(mpas_filename, mpas_source_mesh) /= 0) then
1533             call mprintf(.true.,ERROR, 'Error setting up MPAS mesh %s with scan_input_open', s1=mpas_filename)
1534          end if
1535       end if
1536
1537       !
1538       ! If we have not already defined the WRF grid, do that now...
1539       !
1540       if (.not. wrf_target_mesh_m % valid) then
1541          allocate(xlat_rad(size(xlat,1), size(xlat,2)))
1542          allocate(xlon_rad(size(xlat,1), size(xlat,2)))
1543          xlat_rad(:,:) = xlat(:,:) * deg2rad
1544          xlon_rad(:,:) = xlon(:,:) * deg2rad
1545          call mprintf(.true.,LOGFILE,'Need to set up WRF target mass-grid')
1546          if (target_mesh_setup(wrf_target_mesh_m, lat2d=xlat_rad, lon2d=xlon_rad) /= 0) then
1547             call mprintf(.true.,ERROR, 'Error setting up WRF target grid')
1548          end if
1549
1550          call mprintf(.true.,LOGFILE,'Also computing remapping weights...')
1551          if (remap_info_setup(mpas_source_mesh, wrf_target_mesh_m, remap_info_m) /= 0) then
1552             call mprintf(.true.,ERROR, 'Error computing remapping weights from MPAS to WRF grid')
1553          end if
1554       else
1555          call mprintf(.true.,LOGFILE,'Already set up WRF target mass-grid')
1556       end if
1557
1558       if (.not. wrf_target_mesh_u % valid) then
1559          allocate(xlat_u_rad(size(xlat_u,1), size(xlat_u,2)))
1560          allocate(xlon_u_rad(size(xlat_u,1), size(xlat_u,2)))
1561          xlat_u_rad(:,:) = xlat_u(:,:) * deg2rad
1562          xlon_u_rad(:,:) = xlon_u(:,:) * deg2rad
1563          call mprintf(.true.,LOGFILE,'Need to set up WRF target U-grid')
1564          if (target_mesh_setup(wrf_target_mesh_u, lat2d=xlat_u_rad, lon2d=xlon_u_rad) /= 0) then
1565             call mprintf(.true.,ERROR, 'Error setting up WRF target grid')
1566          end if
1567
1568          call mprintf(.true.,LOGFILE,'Also computing remapping weights...')
1569          if (remap_info_setup(mpas_source_mesh, wrf_target_mesh_u, remap_info_u) /= 0) then
1570             call mprintf(.true.,ERROR, 'Error computing remapping weights from MPAS to WRF grid')
1571          end if
1572       else
1573          call mprintf(.true.,LOGFILE,'Already set up WRF target U-grid')
1574       end if
1575
1576       if (.not. wrf_target_mesh_v % valid) then
1577          allocate(xlat_v_rad(size(xlat_v,1), size(xlat_v,2)))
1578          allocate(xlon_v_rad(size(xlat_v,1), size(xlat_v,2)))
1579          xlat_v_rad(:,:) = xlat_v(:,:) * deg2rad
1580          xlon_v_rad(:,:) = xlon_v(:,:) * deg2rad
1581          call mprintf(.true.,LOGFILE,'Need to set up WRF target V-grid')
1582          if (target_mesh_setup(wrf_target_mesh_v, lat2d=xlat_v_rad, lon2d=xlon_v_rad) /= 0) then
1583             call mprintf(.true.,ERROR, 'Error setting up WRF target grid')
1584          end if
1585
1586          call mprintf(.true.,LOGFILE,'Also computing remapping weights...')
1587          if (remap_info_setup(mpas_source_mesh, wrf_target_mesh_v, remap_info_v) /= 0) then
1588             call mprintf(.true.,ERROR, 'Error computing remapping weights from MPAS to WRF grid')
1589          end if
1590       else
1591          call mprintf(.true.,LOGFILE,'Already set up WRF target V-grid')
1592       end if
1593
1594
1595       if (scan_input_open(mpas_filename, mpas_handle, nRecords) /= 0) then
1596          call mprintf(.true.,ERROR, 'Error opening %s with scan_input_open', s1=mpas_filename)
1597       end if
1598
1599
1600       ! Initialize fg_input structure to store the field
1601       field_to_store%header%version = 1
1602       field_to_store%header%date = '?'
1603       if (do_const_processing) then
1604          field_to_store%header%time_dependent = .false.
1605          field_to_store%header%constant_field = .true.
1606       else
1607          field_to_store%header%time_dependent = .true.
1608          field_to_store%header%constant_field = .false.
1609       end if
1610       field_to_store%header%forecast_hour = 0.0
1611       field_to_store%header%fg_source = 'MPAS'
1612       field_to_store%header%field = ' '
1613       field_to_store%header%field(1:9) = '?'
1614       field_to_store%header%units = ' '
1615       field_to_store%header%units(1:25) = '?'
1616       field_to_store%header%description = ' '
1617       field_to_store%header%description(1:46) = '?'
1618       field_to_store%header%vertical_coord = 'z_dim_name'
1619       field_to_store%header%vertical_level = 0
1620       field_to_store%header%sr_x = 1
1621       field_to_store%header%sr_y = 1
1622       field_to_store%header%array_order = 'XY '
1623       field_to_store%header%is_wind_grid_rel = .false.
1624       field_to_store%header%array_has_missing_values = .false.
1625       nullify(field_to_store%r_arr)
1626       nullify(field_to_store%valid_mask)
1627       nullify(field_to_store%modified_mask)
1628
1629       ! If we should not output this field, just list it as a mask field
1630 !???      if (output_this_field(idx)) then
1631          field_to_store%header%mask_field = .false.
1632 !???      else
1633 !???         field%header%mask_field = .true.
1634 !???      end if
1635
1636
1637       do while (scan_input_next_field(mpas_handle, mpas_field) == 0)
1638
1639          if (can_remap_field(mpas_field)) then
1640
1641             ! Here, rename a few MPAS fields that would be difficult to treat
1642             ! with METGRID.TBL options; principally, these are surface fields
1643             ! that have different names from their upper-air counterparts.
1644             if (trim(mpas_field % name) == 'u10') then
1645                mpas_field % name = 'uReconstructZonal'
1646             else if (trim(mpas_field % name) == 'v10') then
1647                mpas_field % name = 'uReconstructMeridional'
1648             else if (trim(mpas_field % name) == 'q2') then
1649                mpas_field % name = 'qv'
1650             else if (trim(mpas_field % name) == 't2m') then
1651                mpas_field % name = 'theta'
1652             end if
1653
1654             ! Mark this MPAS field as "gotten" for any later checks
1655             ! on mandatory fields
1656             idx = mpas_name_to_idx(trim(mpas_field % name))
1657             if (idx > 0) then
1658                got_this_field(idx) = .true.
1659                if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
1660                   output_flags(idx) = flag_in_output(idx)
1661                end if
1662             else
1663                istat = scan_input_free_field(mpas_field)
1664                cycle
1665             end if
1666
1667             istat = scan_input_read_field(mpas_field, frame=1)
1668
1669             field_to_store%map%stagger = mpas_output_stagger(mpas_field % name)
1670             if (field_to_store%map%stagger == M) then
1671                field_to_store%header%dim1(1) = we_mem_s
1672                field_to_store%header%dim1(2) = we_mem_e
1673                field_to_store%header%dim2(1) = sn_mem_s
1674                field_to_store%header%dim2(2) = sn_mem_e
1675                if (idx > 0) then
1676                   if (masked(idx) == MASKED_WATER) then
1677                      istat = remap_field(remap_info_m, mpas_field, wrf_field, masked=.true.)
1678                   else
1679                      istat = remap_field(remap_info_m, mpas_field, wrf_field, masked=.false.)
1680                   end if
1681                else
1682                   istat = remap_field(remap_info_m, mpas_field, wrf_field, masked=.false.)
1683                end if
1684             else if (field_to_store%map%stagger == U) then
1685                field_to_store%header%dim1(1) = we_mem_stag_s
1686                field_to_store%header%dim1(2) = we_mem_stag_e
1687                field_to_store%header%dim2(1) = sn_mem_s
1688                field_to_store%header%dim2(2) = sn_mem_e
1689                istat = remap_field(remap_info_u, mpas_field, wrf_field)
1690             else if (field_to_store%map%stagger == V) then
1691                field_to_store%header%dim1(1) = we_mem_s
1692                field_to_store%header%dim1(2) = we_mem_e
1693                field_to_store%header%dim2(1) = sn_mem_stag_s
1694                field_to_store%header%dim2(2) = sn_mem_stag_e
1695                istat = remap_field(remap_info_v, mpas_field, wrf_field)
1696             else
1697                call mprintf(.true.,ERROR, 'Cannot handle requested output stagger %i for MPAS field %s ...', &
1698                             i1=field_to_store%map%stagger, s1=trim(mpas_field % name))
1699             end if
1700
1701             if (wrf_field % ndims == 3 .and. trim(wrf_field % dimnames(1)) == 'nVertLevels') then   ! 3-d MPAS atmosphere field
1702                field_to_store % header % field = mpas_to_intermediate_name(mpas_field % name)
1703
1704                ! If no match in the METGRID.TBL was found for this MPAS field, just use the MPAS name
1705                if (len_trim(field_to_store % header % field) == 0) then
1706                   field_to_store % header % field = trim(mpas_field % name)
1707                end if
1708
1709                field_to_store % header % vertical_coord = 'num_mpas_levels'
1710                do k=1,wrf_field % dimlens(1)
1711                   allocate(field_to_store % r_arr(wrf_field % dimlens(2), wrf_field % dimlens(3)))
1712                   allocate(field_to_store % valid_mask)
1713                   call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(2), wrf_field % dimlens(3))
1714                   do j=1,wrf_field % dimlens(3)
1715                   do i=1,wrf_field % dimlens(2)
1716                      call bitarray_set(field_to_store % valid_mask, i, j)
1717                   end do
1718                   end do
1719                   field_to_store % header % vertical_level = k
1720
1721                   if (wrf_field % xtype == FIELD_TYPE_REAL) then
1722                      field_to_store % r_arr(:,:) = wrf_field % array3r(k,:,:)
1723                   else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then
1724                      field_to_store % r_arr(:,:) = wrf_field % array3d(k,:,:)
1725                   end if
1726
1727                   ! The u_field and v_field are used later by calling code to
1728                   ! determine which fields represent the x- and y-components of
1729                   ! horizonal wind velocity for the purposes of wind rotation
1730                   if (trim(mpas_field % name) == 'uReconstructZonal') then
1731                      call dup(field_to_store, u_field)
1732                   end if
1733                   if (trim(mpas_field % name) == 'uReconstructMeridional') then
1734                      call dup(field_to_store, v_field)
1735                   end if
1736
1737                   call storage_put_field(field_to_store)
1738                end do
1739
1740             else if (wrf_field % ndims == 3 .and. trim(wrf_field % dimnames(1)) == 'nVertLevelsP1') then   ! 3-d MPAS atmosphere field
1741                field_to_store % header % field = mpas_to_intermediate_name(mpas_field % name)
1742
1743                ! If no match in the METGRID.TBL was found for this MPAS field, just use the MPAS name
1744                if (len_trim(field_to_store % header % field) == 0) then
1745                   field_to_store % header % field = trim(mpas_field % name)
1746                end if
1747
1748                ! Handle surface level
1749                field_to_store % header % vertical_coord = 'none'
1750                allocate(field_to_store % r_arr(wrf_field % dimlens(2), wrf_field % dimlens(3)))
1751                allocate(field_to_store % valid_mask)
1752                call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(2), wrf_field % dimlens(3))
1753                do j=1,wrf_field % dimlens(3)
1754                do i=1,wrf_field % dimlens(2)
1755                   call bitarray_set(field_to_store % valid_mask, i, j)
1756                end do
1757                end do
1758                field_to_store % header % vertical_level = 200100.0
1759
1760                if (wrf_field % xtype == FIELD_TYPE_REAL) then
1761                   field_to_store % r_arr(:,:) = wrf_field % array3r(1,:,:)
1762                else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then
1763                   field_to_store % r_arr(:,:) = wrf_field % array3d(1,:,:)
1764                end if
1765
1766                call storage_put_field(field_to_store)
1767
1768                ! Handle all other layers
1769                field_to_store % header % vertical_coord = 'num_mpas_levels'
1770                do k=1,wrf_field % dimlens(1) - 1
1771                   allocate(field_to_store % r_arr(wrf_field % dimlens(2), wrf_field % dimlens(3)))
1772                   allocate(field_to_store % valid_mask)
1773                   call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(2), wrf_field % dimlens(3))
1774                   do j=1,wrf_field % dimlens(3)
1775                   do i=1,wrf_field % dimlens(2)
1776                      call bitarray_set(field_to_store % valid_mask, i, j)
1777                   end do
1778                   end do
1779                   field_to_store % header % vertical_level = k
1780
1781                   ! Average to layer midpoint
1782                   if (wrf_field % xtype == FIELD_TYPE_REAL) then
1783                      field_to_store % r_arr(:,:) = 0.5 * (wrf_field % array3r(k,:,:) + wrf_field % array3r(k+1,:,:))
1784                   else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then
1785                      field_to_store % r_arr(:,:) = 0.5 * (wrf_field % array3d(k,:,:) + wrf_field % array3d(k+1,:,:))
1786                   end if
1787
1788                   call storage_put_field(field_to_store)
1789                end do
1790
1791                ! Special case: zgrid field also provides SOILHGT
1792                if (trim(mpas_field % name) == 'zgrid') then
1793                   field_to_store % header % field = 'SOILHGT'
1794
1795                   allocate(field_to_store % r_arr(wrf_field % dimlens(2), wrf_field % dimlens(3)))
1796                   allocate(field_to_store % valid_mask)
1797                   call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(2), wrf_field % dimlens(3))
1798                   do j=1,wrf_field % dimlens(3)
1799                   do i=1,wrf_field % dimlens(2)
1800                      call bitarray_set(field_to_store % valid_mask, i, j)
1801                   end do
1802                   end do
1803                   field_to_store % header % vertical_level = 200100.0
1804
1805                   if (wrf_field % xtype == FIELD_TYPE_REAL) then
1806                      field_to_store % r_arr(:,:) = wrf_field % array3r(1,:,:)
1807                   else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then
1808                      field_to_store % r_arr(:,:) = wrf_field % array3d(1,:,:)
1809                   end if
1810
1811                   call storage_put_field(field_to_store)
1812
1813                   do idx=1,num_entries
1814                      if (trim(fieldname(idx)) == 'SOILHGT') then
1815                         got_this_field(idx) = .true.
1816                         if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
1817                            output_flags(idx) = flag_in_output(idx)
1818                         end if
1819                         exit
1820                      end if
1821                   end do
1822                end if
1823
1824             else if (wrf_field % ndims == 3 .and. trim(wrf_field % dimnames(1)) == 'nSoilLevels') then   ! 3-d MPAS soil field
1825
1826                field_to_store % header % vertical_coord = 'none'
1827                if (trim(mpas_field % name) == 'tslb') then
1828                   do k=1,wrf_field % dimlens(1)
1829                      if (k == 1) then
1830                         field_to_store % header % field = 'ST000010'
1831                      else if (k == 2) then
1832                         field_to_store % header % field = 'ST010040'
1833                      else if (k == 3) then
1834                         field_to_store % header % field = 'ST040100'
1835                      else if (k == 4) then
1836                         field_to_store % header % field = 'ST100200'
1837                      else
1838                         call mprintf(.true.,ERROR, 'Too many soil layers in MPAS soil field %s ...', s1=trim(mpas_field % name))
1839                      end if
1840                      allocate(field_to_store % r_arr(wrf_field % dimlens(2), wrf_field % dimlens(3)))
1841                      allocate(field_to_store % valid_mask)
1842                      call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(2), wrf_field % dimlens(3))
1843                      do j=1,wrf_field % dimlens(3)
1844                      do i=1,wrf_field % dimlens(2)
1845                         call bitarray_set(field_to_store % valid_mask, i, j)
1846                      end do
1847                      end do
1848                      field_to_store % header % vertical_level = 200100.0
1849
1850                      if (wrf_field % xtype == FIELD_TYPE_REAL) then
1851                         field_to_store % r_arr(:,:) = wrf_field % array3r(k,:,:)
1852                      else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then
1853                         field_to_store % r_arr(:,:) = wrf_field % array3d(k,:,:)
1854                      end if
1855
1856                      if (idx > 0) then
1857                         if (masked(idx) == MASKED_WATER) then
1858                            where (landmask(:,:) == 0) field_to_store % r_arr(:,:) = fill_missing(idx)
1859                         end if
1860                      end if
1861
1862                      call storage_put_field(field_to_store)
1863                   end do
1864                else if (trim(mpas_field % name) == 'smois') then
1865                   do k=1,wrf_field % dimlens(1)
1866                      if (k == 1) then
1867                         field_to_store % header % field = 'SM000010'
1868                      else if (k == 2) then
1869                         field_to_store % header % field = 'SM010040'
1870                      else if (k == 3) then
1871                         field_to_store % header % field = 'SM040100'
1872                      else if (k == 4) then
1873                         field_to_store % header % field = 'SM100200'
1874                      else
1875                         call mprintf(.true.,ERROR, 'Too many soil layers in MPAS soil field %s ...', s1=trim(mpas_field % name))
1876                      end if
1877                      allocate(field_to_store % r_arr(wrf_field % dimlens(2), wrf_field % dimlens(3)))
1878                      allocate(field_to_store % valid_mask)
1879                      call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(2), wrf_field % dimlens(3))
1880                      do j=1,wrf_field % dimlens(3)
1881                      do i=1,wrf_field % dimlens(2)
1882                         call bitarray_set(field_to_store % valid_mask, i, j)
1883                      end do
1884                      end do
1885                      field_to_store % header % vertical_level = 200100.0
1886
1887                      if (wrf_field % xtype == FIELD_TYPE_REAL) then
1888                         field_to_store % r_arr(:,:) = wrf_field % array3r(k,:,:)
1889                      else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then
1890                         field_to_store % r_arr(:,:) = wrf_field % array3d(k,:,:)
1891                      end if
1892
1893                      if (idx > 0) then
1894                         if (masked(idx) == MASKED_WATER) then
1895                            where (landmask(:,:) == 0) field_to_store % r_arr(:,:) = fill_missing(idx)
1896                         end if
1897                      end if
1898
1899                      call storage_put_field(field_to_store)
1900                   end do
1901                else
1902                   call mprintf(.true.,WARN, 'Skipping unknown MPAS soil field %s ...', s1=trim(mpas_field % name))
1903                end if
1904
1905             else if (wrf_field % ndims == 2) then   ! 2-d MPAS field
1906                field_to_store % header % field = mpas_to_intermediate_name(mpas_field % name)
1907
1908                ! If no match in the METGRID.TBL was found for this MPAS field, just use the MPAS name
1909                if (len_trim(field_to_store % header % field) == 0) then
1910                   field_to_store % header % field = trim(mpas_field % name)
1911                end if
1912
1913                field_to_store % header % vertical_coord = 'none'
1914                allocate(field_to_store % r_arr(wrf_field % dimlens(1), wrf_field % dimlens(2)))
1915                allocate(field_to_store % valid_mask)
1916                call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(1), wrf_field % dimlens(2))
1917                do j=1,wrf_field % dimlens(2)
1918                do i=1,wrf_field % dimlens(1)
1919                   call bitarray_set(field_to_store % valid_mask, i, j)
1920                end do
1921                end do
1922                field_to_store % header % vertical_level = 200100.0
1923
1924                if (wrf_field % xtype == FIELD_TYPE_REAL) then
1925                   field_to_store % r_arr(:,:) = wrf_field % array2r(:,:)
1926                else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then
1927                   field_to_store % r_arr(:,:) = wrf_field % array2d(:,:)
1928                end if
1929
1930                if (idx > 0) then
1931                   if (masked(idx) == MASKED_WATER) then
1932                      where (landmask(:,:) == 0) field_to_store % r_arr(:,:) = fill_missing(idx)
1933                   end if
1934                end if
1935
1936                call storage_put_field(field_to_store)
1937             end if
1938
1939             istat = free_target_field(wrf_field)
1940          end if
1941
1942          istat = scan_input_free_field(mpas_field)
1943       end do
1944
1945       if (scan_input_close(mpas_handle) /= 0) then
1946          call mprintf(.true.,ERROR, 'Error closing %s with scan_input_close', s1=mpas_filename)
1947       end if
1948
1949    end subroutine process_mpas_fields
1950
1951
1952    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1953    ! Name: derive_mpas_fields
1954    !
1955    ! Purpose:
1956    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1957    subroutine derive_mpas_fields()
1958
1959       use bitarray_module
1960       use gridinfo_module
1961       use interp_module
1962       use interp_option_module
1963       use read_met_module
1964       use llxy_module
1965       use misc_definitions_module
1966       use module_debug
1967       use output_module
1968       use parallel_module
1969       use rotate_winds_module
1970       use storage_module
1971       use scan_input
1972       use mpas_mesh
1973       use constants_module
1974
1975       implicit none
1976
1977       integer :: k
1978       integer :: istatus
1979       integer, pointer, dimension(:) :: theta_levels => null()
1980       integer, pointer, dimension(:) :: pressure_levels => null()
1981       real, pointer, dimension(:,:) :: exner => null()
1982       type (fg_input) :: theta_field
1983       type (fg_input) :: pressure_field
1984
1985       !
1986       ! Derive TT from theta and pressure
1987       !
1988       theta_field%header%time_dependent = .true.
1989       theta_field%header%constant_field = .false.
1990       theta_field%header%field = 'TT'
1991       theta_field%header%vertical_level = 0
1992       nullify(theta_field%r_arr)
1993       nullify(theta_field%valid_mask)
1994       nullify(theta_field%modified_mask)
1995
1996       pressure_field%header%time_dependent = .true.
1997       pressure_field%header%constant_field = .false.
1998       pressure_field%header%field = 'PRESSURE'
1999       pressure_field%header%vertical_level = 0
2000       nullify(pressure_field%r_arr)
2001       nullify(pressure_field%valid_mask)
2002       nullify(pressure_field%modified_mask)
2003
2004       call storage_get_levels(theta_field, theta_levels)
2005       call storage_get_levels(pressure_field, pressure_levels)
2006
2007       if (associated(theta_levels) .and. associated(pressure_levels)) then
2008 !         call mprintf(.true.,LOGFILE, 'Computing MPAS TT field from theta and pressure...')
2009
2010          if (size(theta_levels) == size(pressure_levels)) then
2011             do k = 1, size(theta_levels)
2012                theta_field % header % vertical_level = theta_levels(k)
2013                call storage_get_field(theta_field, istatus)
2014                if (trim(theta_field % header % fg_source) /= 'MPAS') then
2015                   cycle
2016                end if
2017                if (istatus /= 0) then
2018                   call mprintf(.true.,ERROR, 'Could not get MPAS theta field at level %i', i1=theta_levels(k))
2019                   return
2020                end if
2021
2022                pressure_field % header % vertical_level = pressure_levels(k)
2023                call storage_get_field(pressure_field, istatus)
2024                if (trim(pressure_field % header % fg_source) /= 'MPAS') then
2025                   cycle
2026                end if
2027                if (istatus /= 0) then
2028                   call mprintf(.true.,ERROR, 'Could not get MPAS pressure field at level %i', i1=theta_levels(k))
2029                   return
2030                end if
2031
2032                ! Compute temperature
2033                call mprintf(.true.,LOGFILE, 'Computing TT at level %i for MPAS dataset', i1=theta_levels(k))
2034                if (.not. associated(exner)) then
2035                   allocate(exner(size(theta_field % r_arr, 1), size(theta_field % r_arr, 2)))
2036                end if
2037                exner(:,:) = (pressure_field % r_arr(:,:) / P0)**(RD/CP)
2038                theta_field % r_arr(:,:) = theta_field % r_arr(:,:) * exner(:,:)
2039
2040                call storage_put_field(theta_field)
2041             end do
2042 !         else
2043 !            call mprintf(.true.,ERROR, 'The MPAS theta and pressure fields do not have the same number of levels!')
2044          end if
2045
2046          deallocate(theta_levels)
2047          deallocate(pressure_levels)
2048          if (associated(exner)) then
2049             deallocate(exner)
2050          end if
2051       end if
2052
2053    end subroutine derive_mpas_fields
2054
2055
2056    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2057    ! Name: get_interp_masks
2058    !
2059    ! Purpose:
2060    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2061    subroutine get_interp_masks(fg_prefix, is_constants, fg_date)
2062
2063       use interp_option_module
2064       use read_met_module
2065       use storage_module
2066
2067       implicit none
2068
2069       ! Arguments
2070       logical, intent(in) :: is_constants
2071       character (len=*), intent(in) :: fg_prefix, fg_date
2072
2073 ! BUG: Move this constant to misc_definitions_module?
2074 integer, parameter :: BDR_WIDTH = 3
2075
2076       ! Local variables
2077       integer :: i, istatus, idx, idxt
2078       type (fg_input) :: mask_field
2079       type (met_data) :: fg_data
2080
2081       istatus = 0
2082
2083       call read_met_init(fg_prefix, is_constants, fg_date, istatus)
2084
2085       do while (istatus == 0)
2086
2087          call read_next_met_field(fg_data, istatus)
2088
2089          if (istatus == 0) then
2090
2091             ! Find out which METGRID.TBL entry goes with this field
2092             idxt = num_entries + 1
2093             do idx=1,num_entries
2094                if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
2095                    (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
2096
2097                   if (index(fg_prefix,trim(from_input(idx))) /= 0 .or. &
2098                      (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
2099                      idxt = idx
2100                   end if
2101
2102                end if
2103             end do
2104             idx = idxt
2105             if (idx > num_entries) idx = num_entries ! The last entry is a default
2106
2107             ! Do we need to rename this field?
2108             if (output_name(idx) /= ' ') then
2109                fg_data%field = output_name(idx)(1:9)
2110
2111                idxt = num_entries + 1
2112                do idx=1,num_entries
2113                   if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
2114                       (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
2115
2116                      if (index(fg_prefix,trim(from_input(idx))) /= 0 .or. &
2117                         (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
2118                         idxt = idx
2119                      end if
2120
2121                   end if
2122                end do
2123                idx = idxt
2124                if (idx > num_entries) idx = num_entries ! The last entry is a default
2125             end if
2126
2127             do i=1,num_entries
2128                if (interp_mask(i) /= ' ' .and. (trim(interp_mask(i)) == trim(fg_data%field))) then
2129
2130                   mask_field%header%version = 1
2131                   mask_field%header%date = ' '
2132                   mask_field%header%date = fg_date
2133                   if (is_constants) then
2134                      mask_field%header%time_dependent = .false.
2135                      mask_field%header%constant_field = .true.
2136                   else
2137                      mask_field%header%time_dependent = .true.
2138                      mask_field%header%constant_field = .false.
2139                   end if
2140                   mask_field%header%mask_field = .true.
2141                   mask_field%header%forecast_hour = 0.
2142                   mask_field%header%fg_source = 'degribbed met data'
2143                   mask_field%header%field = trim(fg_data%field)//'.mask'
2144                   mask_field%header%units = '-'
2145                   mask_field%header%description = '-'
2146                   mask_field%header%vertical_coord = 'none'
2147                   mask_field%header%vertical_level = 1
2148                   mask_field%header%sr_x = 1
2149                   mask_field%header%sr_y = 1
2150                   mask_field%header%array_order = 'XY'
2151                   mask_field%header%dim1(1) = 1
2152                   mask_field%header%dim1(2) = fg_data%nx
2153                   mask_field%header%dim2(1) = 1
2154                   mask_field%header%dim2(2) = fg_data%ny
2155                   mask_field%header%is_wind_grid_rel = .true.
2156                   mask_field%header%array_has_missing_values = .false.
2157                   mask_field%map%stagger = M
2158
2159                   ! Do a simple check to see whether this is a global lat/lon dataset
2160                   ! Note that we do not currently support regional Gaussian grids
2161                   if ((fg_data%iproj == PROJ_LATLON .and. abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) &
2162                        .or. (fg_data%iproj == PROJ_GAUSS)) then
2163                      allocate(mask_field%r_arr(1-BDR_WIDTH:fg_data%nx+BDR_WIDTH,1:fg_data%ny))
2164
2165                      mask_field%r_arr(1:fg_data%nx,                      1:fg_data%ny) = &
2166                          fg_data%slab(1:fg_data%nx,              1:fg_data%ny)
2167
2168                      mask_field%r_arr(1-BDR_WIDTH:0,                     1:fg_data%ny) = &
2169                          fg_data%slab(fg_data%nx-BDR_WIDTH+1:fg_data%nx, 1:fg_data%ny)
2170
2171                      mask_field%r_arr(fg_data%nx+1:fg_data%nx+BDR_WIDTH, 1:fg_data%ny) = &
2172                          fg_data%slab(1:BDR_WIDTH,       1:fg_data%ny)
2173
2174                   else
2175                      allocate(mask_field%r_arr(1:fg_data%nx,1:fg_data%ny))
2176                      mask_field%r_arr = fg_data%slab
2177                   end if
2178
2179                   nullify(mask_field%valid_mask)
2180                   nullify(mask_field%modified_mask)
2181
2182                   call storage_put_field(mask_field)
2183
2184                   exit
2185
2186                end if
2187             end do
2188
2189             if (associated(fg_data%slab)) deallocate(fg_data%slab)
2190
2191          end if
2192
2193       end do
2194
2195       call read_met_close()
2196
2197    end subroutine get_interp_masks
2198
2199
2200    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2201    ! Name: interp_met_field
2202    !
2203    ! Purpose:
2204    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2205    subroutine interp_met_field(input_name, short_fieldnm, ifieldstagger, istagger, &
2206                                field, xlat, xlon, sm1, em1, sm2, em2, &
2207                                sd1, ed1, sd2, ed2, &
2208                                slab, minx, maxx, miny, maxy, bdr, do_gcell_interp, &
2209                                new_pts, process_bdy_width, landmask)
2210
2211       use bitarray_module
2212       use interp_module
2213       use interp_option_module
2214       use llxy_module
2215       use misc_definitions_module
2216       use storage_module
2217
2218       implicit none
2219
2220       ! Arguments
2221       integer, intent(in) :: ifieldstagger, istagger, &
2222                              sm1, em1, sm2, em2, &
2223                              sd1, ed1, sd2, ed2, &
2224                              minx, maxx, miny, maxy, bdr, &
2225                              process_bdy_width
2226       real, dimension(minx:maxx,miny:maxy), intent(in) :: slab
2227       real, dimension(sm1:em1,sm2:em2), intent(in) :: xlat, xlon
2228       real, dimension(sm1:em1,sm2:em2), intent(in), optional :: landmask
2229       logical, intent(in) :: do_gcell_interp
2230       character (len=9), intent(in) :: short_fieldnm
2231       character (len=MAX_FILENAME_LEN), intent(in) :: input_name
2232       type (fg_input), intent(inout) :: field
2233       type (bitarray), intent(inout) :: new_pts
2234
2235       ! Local variables
2236       integer :: i, j, idx, idxt, orig_selected_proj, interp_mask_status, &
2237                  interp_land_mask_status, interp_water_mask_status, process_width
2238       integer, pointer, dimension(:) :: interp_array, interp_opts
2239       real :: rx, ry, temp
2240       real, pointer, dimension(:,:) :: data_count
2241       type (fg_input) :: mask_field, mask_water_field, mask_land_field
2242       !BPR BEGIN
2243       real, dimension(sm1:em1,sm2:em2) :: r_arr_cur_source
2244       !BPR END
2245
2246       ! CWH Initialize local pointer variables
2247       nullify(interp_array)
2248       nullify(interp_opts)
2249       nullify(data_count)
2250
2251       ! Find index into fieldname, interp_method, masked, and fill_missing
2252       !   of the current field
2253       idxt = num_entries + 1
2254       do idx=1,num_entries
2255          if ((index(fieldname(idx), trim(short_fieldnm)) /= 0) .and. &
2256              (len_trim(fieldname(idx)) == len_trim(short_fieldnm))) then
2257             if (index(input_name,trim(from_input(idx))) /= 0 .or. &
2258                (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
2259                idxt = idx
2260             end if
2261          end if
2262       end do
2263       idx = idxt
2264       if (idx > num_entries) then
2265          call mprintf(.true.,WARN,'Entry in METGRID.TBL not found for field %s. '// &
2266                       'Default options will be used for this field!', s1=short_fieldnm)
2267          idx = num_entries ! The last entry is a default
2268       end if
2269
2270       if (process_bdy_width == 0) then
2271          process_width = max(ed1-sd1+1, ed2-sd2+1)
2272       else
2273          process_width = process_bdy_width
2274         ! Add two extra rows/cols to accommodate staggered fields: one extra row/col for
2275         !    averaging to mass points in real, and one beyond that for averaging during
2276         !    wind rotation
2277         if (ifieldstagger /= M) process_width = process_width + 2
2278       end if
2279
2280       field%header%dim1(1) = sm1
2281       field%header%dim1(2) = em1
2282       field%header%dim2(1) = sm2
2283       field%header%dim2(2) = em2
2284       field%map%stagger = ifieldstagger
2285       if (.not. associated(field%r_arr)) then
2286          allocate(field%r_arr(sm1:em1,sm2:em2))
2287       end if
2288
2289       interp_mask_status = 1
2290       interp_land_mask_status = 1
2291       interp_water_mask_status = 1
2292
2293       if (interp_mask(idx) /= ' ') then
2294          mask_field%header%version = 1
2295          mask_field%header%forecast_hour = 0.
2296          mask_field%header%field = trim(interp_mask(idx))//'.mask'
2297          mask_field%header%vertical_coord = 'none'
2298          mask_field%header%vertical_level = 1
2299
2300          call storage_get_field(mask_field, interp_mask_status)
2301
2302       end if
2303       if (interp_land_mask(idx) /= ' ') then
2304          mask_land_field%header%version = 1
2305          mask_land_field%header%forecast_hour = 0.
2306          mask_land_field%header%field = trim(interp_land_mask(idx))//'.mask'
2307          mask_land_field%header%vertical_coord = 'none'
2308          mask_land_field%header%vertical_level = 1
2309
2310          call storage_get_field(mask_land_field, interp_land_mask_status)
2311
2312       end if
2313       if (interp_water_mask(idx) /= ' ') then
2314          mask_water_field%header%version = 1
2315          mask_water_field%header%forecast_hour = 0.
2316          mask_water_field%header%field = trim(interp_water_mask(idx))//'.mask'
2317          mask_water_field%header%vertical_coord = 'none'
2318          mask_water_field%header%vertical_level = 1
2319
2320          call storage_get_field(mask_water_field, interp_water_mask_status)
2321
2322       end if
2323
2324       interp_array => interp_array_from_string(interp_method(idx))
2325       interp_opts => interp_options_from_string(interp_method(idx))
2326
2327
2328       !
2329       ! Interpolate using average_gcell interpolation method
2330       !
2331       if (do_gcell_interp) then
2332          !BPR BEGIN
2333          !If a lower priority source of the current field has already been read
2334          !in, the results of that input are currently in field%r_arr
2335          !Pass COPY of field%r_arr into accum_continous because in accum_continuous
2336          !it will set the input variable to zero over points covered by the
2337          !current source and then determine the appropriate value to place at
2338          !that point.  This will overwrite data already put in field%r_arr by
2339          !lower priority sources.
2340          !This is problematic if the current source results in missing values
2341          !over parts of the area covered by the current source where a lower
2342          !priority source has already provided a value. In this case, if one
2343          !passes in field%r_arr, it will overwrite the value provided by the
2344          !lower priority source with zero.
2345          r_arr_cur_source = field%r_arr
2346          !BPR END
2347          allocate(data_count(sm1:em1,sm2:em2))
2348          data_count = 0.
2349
2350          if (interp_mask_status == 0) then
2351             call accum_continuous(slab, &
2352                          minx, maxx, miny, maxy, 1, 1, bdr, &
2353                          !BPR BEGIN
2354                          !Pass COPY of field%r_arr instead of field%r_arr itself
2355                          !field%r_arr, data_count, &
2356                          r_arr_cur_source, data_count, &
2357                          !BPR END
2358                          sm1, em1, sm2, em2, 1, 1, &
2359                          istagger, &
2360                          new_pts, missing_value(idx), interp_mask_val(idx), interp_mask_relational(idx), mask_field%r_arr)
2361          else
2362             call accum_continuous(slab, &
2363                          minx, maxx, miny, maxy, 1, 1, bdr, &
2364                          !BPR BEGIN
2365                          !Pass COPY of field%r_arr instead of field%r_arr itself
2366                          !field%r_arr, data_count, &
2367                          r_arr_cur_source, data_count, &
2368                          !BPR END
2369                          sm1, em1, sm2, em2, 1, 1, &
2370                          istagger, &
2371                          new_pts, missing_value(idx), -1.) ! The -1 is the maskval, but since we
2372                                                            !   we do not give an optional mask, no
2373                                                            !   no need to worry about -1s in data
2374          end if
2375
2376          orig_selected_proj = iget_selected_domain()
2377          call select_domain(SOURCE_PROJ)
2378          do j=sm2,em2
2379             do i=sm1,em1
2380
2381                if (abs(i - sd1) >= process_width .and. (abs(i - ed1) >= process_width) .and. &
2382                    abs(j - sd2) >= process_width .and. (abs(j - ed2) >= process_width)) then
2383                   field%r_arr(i,j) = fill_missing(idx)
2384                   call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2385                   cycle
2386                end if
2387
2388                if (present(landmask)) then
2389
2390                   if (landmask(i,j) /= masked(idx)) then
2391                      if (data_count(i,j) > 0.) then
2392                         !BPR BEGIN
2393                         !accum_continuous is now passed a copy of field%r_arr (r_arr_cur_source)
2394                         !instead of field%r_arr itself so that it does not set
2395                         !field%r_arr to zero where the input source is marked as missing
2396                         !field%r_arr(i,j) = field%r_arr(i,j) / data_count(i,j)
2397                         field%r_arr(i,j) = r_arr_cur_source(i,j) / data_count(i,j)
2398                         !BPR END
2399                         call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2400                      else
2401
2402                         if (interp_mask_status == 0) then
2403                            temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2404                                                    minx, maxx, miny, maxy, bdr, missing_value(idx), &
2405                                                    mask_relational=interp_mask_relational(idx), &
2406                                                    mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
2407                         else
2408                            temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2409                                                    minx, maxx, miny, maxy, bdr, missing_value(idx))
2410                         end if
2411
2412                         if (temp /= missing_value(idx)) then
2413                            field%r_arr(i,j) = temp
2414                            call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2415                         end if
2416
2417                      end if
2418                   else
2419                      field%r_arr(i,j) = fill_missing(idx)
2420                      call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2421                   end if
2422
2423                   if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
2424                       .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
2425                      field%r_arr(i,j) = fill_missing(idx)
2426
2427                      ! Assume that if missing fill value is other than default, then user has asked
2428                      !    to fill in any missing values, and we can consider this point to have
2429                      !    received a valid value
2430                      if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2431                   end if
2432
2433                else
2434
2435                   if (data_count(i,j) > 0.) then
2436                      !BPR BEGIN
2437                      !accum_continuous is now passed a copy of field%r_arr (r_arr_cur_source)
2438                      !instead of field%r_arr itself so that it does not set
2439                      !field%r_arr to zero where the input source is marked as missing
2440                      !field%r_arr(i,j) = field%r_arr(i,j) / data_count(i,j)
2441                      field%r_arr(i,j) = r_arr_cur_source(i,j) / data_count(i,j)
2442                      !BPR END
2443                      call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2444                   else
2445
2446                      if (interp_mask_status == 0) then
2447                         temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2448                                                 minx, maxx, miny, maxy, bdr, missing_value(idx), &
2449                                                 mask_relational=interp_mask_relational(idx), &
2450                                                 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
2451                      else
2452                         temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2453                                                 minx, maxx, miny, maxy, bdr, missing_value(idx))
2454                      end if
2455
2456                      if (temp /= missing_value(idx)) then
2457                         field%r_arr(i,j) = temp
2458                         call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2459                      end if
2460
2461                   end if
2462
2463                   if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
2464                       .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
2465                      field%r_arr(i,j) = fill_missing(idx)
2466
2467                      ! Assume that if missing fill value is other than default, then user has asked
2468                      !    to fill in any missing values, and we can consider this point to have
2469                      !    received a valid value
2470                      if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2471                   end if
2472
2473                end if
2474
2475             end do
2476          end do
2477          call select_domain(orig_selected_proj)
2478          deallocate(data_count)
2479
2480       !
2481       ! No average_gcell interpolation method
2482       !
2483       else
2484
2485          orig_selected_proj = iget_selected_domain()
2486          call select_domain(SOURCE_PROJ)
2487          do j=sm2,em2
2488             do i=sm1,em1
2489
2490                if (abs(i - sd1) >= process_width .and. (abs(i - ed1) >= process_width) .and. &
2491                    abs(j - sd2) >= process_width .and. (abs(j - ed2) >= process_width)) then
2492                   field%r_arr(i,j) = fill_missing(idx)
2493                   call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2494                   cycle
2495                end if
2496
2497                if (present(landmask)) then
2498
2499                   if (masked(idx) == MASKED_BOTH) then
2500
2501                      if (landmask(i,j) == 0) then  ! WATER POINT
2502
2503                         if (interp_land_mask_status == 0) then
2504                            temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2505                                                    minx, maxx, miny, maxy, bdr, missing_value(idx), &
2506                                                    mask_relational=interp_land_mask_relational(idx), &
2507                                                    mask_val=interp_land_mask_val(idx), mask_field=mask_land_field%r_arr)
2508                         else
2509                            temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2510                                                    minx, maxx, miny, maxy, bdr, missing_value(idx))
2511                         end if
2512
2513                      else if (landmask(i,j) == 1) then  ! LAND POINT
2514
2515                         if (interp_water_mask_status == 0) then
2516                            temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2517                                                    minx, maxx, miny, maxy, bdr, missing_value(idx), &
2518                                                    mask_relational=interp_water_mask_relational(idx), &
2519                                                    mask_val=interp_water_mask_val(idx), mask_field=mask_water_field%r_arr)
2520                         else
2521                            temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2522                                                    minx, maxx, miny, maxy, bdr, missing_value(idx))
2523                         end if
2524
2525                      end if
2526
2527                   else if (landmask(i,j) /= masked(idx)) then
2528
2529                      if (interp_mask_status == 0) then
2530                         temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2531                                                 minx, maxx, miny, maxy, bdr, missing_value(idx), &
2532                                                 mask_relational=interp_mask_relational(idx), &
2533                                                 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
2534                      else
2535                         temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2536                                                 minx, maxx, miny, maxy, bdr, missing_value(idx))
2537                      end if
2538
2539                   else
2540                      temp = missing_value(idx)
2541                   end if
2542
2543                ! No landmask for this field
2544                else
2545
2546                   if (interp_mask_status == 0) then
2547                      temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2548                                              minx, maxx, miny, maxy, bdr, missing_value(idx), &
2549                                              mask_relational=interp_mask_relational(idx), &
2550                                              mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
2551                   else
2552                      temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2553                                              minx, maxx, miny, maxy, bdr, missing_value(idx))
2554                   end if
2555
2556                end if
2557
2558                if (temp /= missing_value(idx)) then
2559                   field%r_arr(i,j) = temp
2560                   call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2561                else if (present(landmask)) then
2562                   if (landmask(i,j) == masked(idx)) then
2563                      field%r_arr(i,j) = fill_missing(idx)
2564                      call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2565                   end if
2566                end if
2567
2568                if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
2569                    .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
2570                   field%r_arr(i,j) = fill_missing(idx)
2571
2572                   ! Assume that if missing fill value is other than default, then user has asked
2573                   !    to fill in any missing values, and we can consider this point to have
2574                   !    received a valid value
2575                   if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2576                end if
2577
2578             end do
2579          end do
2580          call select_domain(orig_selected_proj)
2581       end if
2582
2583       deallocate(interp_array)
2584       deallocate(interp_opts)
2585
2586    end subroutine interp_met_field
2587
2588
2589    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2590    ! Name: interp_to_latlon
2591    !
2592    ! Purpose:
2593    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2594    function interp_to_latlon(rlat, rlon, istagger, interp_method_list, interp_opt_list, slab, &
2595                              minx, maxx, miny, maxy, bdr, source_missing_value, &
2596                              mask_field, mask_relational, mask_val)
2597
2598       use interp_module
2599       use llxy_module
2600
2601       implicit none
2602
2603       ! Arguments
2604       integer, intent(in) :: minx, maxx, miny, maxy, bdr, istagger
2605       integer, dimension(:), intent(in) :: interp_method_list
2606       integer, dimension(:), intent(in) :: interp_opt_list
2607       real, intent(in) :: rlat, rlon, source_missing_value
2608       real, dimension(minx:maxx,miny:maxy), intent(in) :: slab
2609       real, intent(in), optional :: mask_val
2610       real, dimension(minx:maxx,miny:maxy), intent(in), optional :: mask_field
2611       character(len=1), intent(in), optional :: mask_relational
2612
2613       ! Return value
2614       real :: interp_to_latlon
2615
2616       ! Local variables
2617       real :: rx, ry
2618
2619       interp_to_latlon = source_missing_value
2620
2621       call lltoxy(rlat, rlon, rx, ry, istagger)
2622       if (rx >= minx+bdr-0.5 .and. rx <= maxx-bdr+0.5) then
2623          if (present(mask_field) .and. present(mask_val) .and. present (mask_relational)) then
2624             interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
2625                                                interp_method_list, interp_opt_list, 1, mask_relational, mask_val, mask_field)
2626          else if (present(mask_field) .and. present(mask_val)) then
2627             interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
2628                                                interp_method_list, interp_opt_list, 1, maskval=mask_val, mask_array=mask_field)
2629          else
2630             interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
2631                                                interp_method_list, interp_opt_list, 1)
2632          end if
2633       else
2634          interp_to_latlon = source_missing_value
2635       end if
2636
2637       if (interp_to_latlon == source_missing_value) then
2638
2639          ! Try a lon in the range 0. to 360.; all lons in the xlon
2640          !    array should be in the range -180. to 180.
2641          if (rlon < 0.) then
2642             call lltoxy(rlat, rlon+360., rx, ry, istagger)
2643             if (rx >= minx+bdr-0.5 .and. rx <= maxx-bdr+0.5) then
2644                if (present(mask_field) .and. present(mask_val) .and. present(mask_relational)) then
2645                   interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
2646                                                      1, 1, source_missing_value, &
2647                                                      interp_method_list, interp_opt_list, 1, &
2648                                                      mask_relational, mask_val, mask_field)
2649                else if (present(mask_field) .and. present(mask_val)) then
2650                   interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
2651                                                      1, 1, source_missing_value, &
2652                                                      interp_method_list, interp_opt_list, 1, &
2653                                                      maskval=mask_val, mask_array=mask_field)
2654                else
2655                   interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
2656                                                      1, 1, source_missing_value, &
2657                                                      interp_method_list, interp_opt_list, 1)
2658                end if
2659             else
2660                interp_to_latlon = source_missing_value
2661             end if
2662
2663          end if
2664
2665       end if
2666
2667       return
2668
2669    end function interp_to_latlon
2670
2671
2672    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2673    ! Name: get_bottom_top_dim
2674    !
2675    ! Purpose:
2676    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2677    subroutine get_bottom_top_dim(bottom_top_dim)
2678
2679       use interp_option_module
2680       use list_module
2681       use storage_module
2682
2683       implicit none
2684
2685       ! Arguments
2686       integer, intent(out) :: bottom_top_dim
2687
2688       ! Local variables
2689       integer :: i, j
2690       integer, pointer, dimension(:) :: field_levels
2691       character (len=32) :: z_dim
2692       type (fg_input), pointer, dimension(:) :: headers
2693       type (list) :: temp_levels
2694
2695       !CWH Initialize local pointer variables
2696       nullify(field_levels)
2697       nullify(headers)
2698
2699       ! Initialize a list to store levels that are found for 3-d fields
2700       call list_init(temp_levels)
2701
2702       ! Get a list of all time-dependent fields (given by their headers) from
2703       !   the storage module
2704       call storage_get_td_headers(headers)
2705
2706       !
2707       ! Given headers of all fields, we first build a list of all possible levels
2708       !    for 3-d met fields (excluding sea-level, though).
2709       !
2710       do i=1,size(headers)
2711          call get_z_dim_name(headers(i)%header%field, z_dim)
2712
2713          ! We only want to consider 3-d met fields
2714          if (z_dim(1:18) == 'num_metgrid_levels') then
2715
2716             ! Find out what levels the current field has
2717             call storage_get_levels(headers(i), field_levels)
2718             do j=1,size(field_levels)
2719
2720                ! If this level has not yet been encountered, add it to our list
2721                if (.not. list_search(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))) then
2722                   if (field_levels(j) /= 201300) then
2723                      call list_insert(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))
2724                   end if
2725                end if
2726
2727             end do
2728
2729             deallocate(field_levels)
2730
2731          end if
2732
2733       end do
2734
2735       bottom_top_dim = list_length(temp_levels)
2736
2737       call list_destroy(temp_levels)
2738       deallocate(headers)
2739
2740    end subroutine get_bottom_top_dim
2741
2742
2743    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2744    ! Name: fill_missing_levels
2745    !
2746    ! Purpose:
2747    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2748    subroutine fill_missing_levels(output_flags)
2749
2750       use interp_option_module
2751       use list_module
2752       use module_debug
2753       use module_mergesort
2754       use storage_module
2755
2756       implicit none
2757
2758       ! Arguments
2759       character (len=128), dimension(:), intent(inout) :: output_flags
2760
2761       ! Local variables
2762       integer :: i, ii, j, ix, jx, k, lower, upper, temp, istatus
2763       integer, pointer, dimension(:) :: union_levels, field_levels
2764       real, pointer, dimension(:) :: r_union_levels
2765       character (len=128) :: clevel
2766       type (fg_input) :: lower_field, upper_field, new_field, search_field
2767       type (fg_input), pointer, dimension(:) :: headers, all_headers
2768       type (list) :: temp_levels
2769       type (list_item), pointer, dimension(:) :: keys
2770
2771       ! CWH Initialize local pointer variables
2772       nullify(union_levels)
2773       nullify(field_levels)
2774       nullify(r_union_levels)
2775       nullify(headers)
2776       nullify(all_headers)
2777       nullify(keys)
2778
2779       ! Initialize a list to store levels that are found for 3-d fields
2780       call list_init(temp_levels)
2781
2782       ! Get a list of all fields (given by their headers) from the storage module
2783       call storage_get_td_headers(headers)
2784       call storage_get_all_headers(all_headers)
2785
2786       !
2787       ! Given headers of all fields, we first build a list of all possible levels
2788       !    for 3-d met fields (excluding sea-level, though).
2789       !
2790       do i=1,size(headers)
2791
2792          ! Find out what levels the current field has
2793          call storage_get_levels(headers(i), field_levels)
2794          do j=1,size(field_levels)
2795
2796             ! If this level has not yet been encountered, add it to our list
2797             if (.not. list_search(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))) then
2798                if (field_levels(j) /= 201300) then
2799                   call list_insert(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))
2800                end if
2801             end if
2802
2803          end do
2804
2805          deallocate(field_levels)
2806
2807       end do
2808
2809       if (list_length(temp_levels) > 0) then
2810
2811          !
2812          ! With all possible levels stored in a list, get an array of levels, sorted
2813          !    in decreasing order
2814          !
2815          i = 0
2816          allocate(union_levels(list_length(temp_levels)))
2817          do while (list_length(temp_levels) > 0)
2818             i = i + 1
2819             call list_get_first_item(temp_levels, ikey=union_levels(i), ivalue=temp)
2820          end do
2821          call mergesort(union_levels, 1, size(union_levels))
2822
2823          allocate(r_union_levels(size(union_levels)))
2824          do i=1,size(union_levels)
2825             r_union_levels(i) = real(union_levels(i))
2826          end do
2827
2828          !
2829          ! With a sorted, complete list of levels, we need
2830          !    to go back and fill in missing levels for each 3-d field
2831          !
2832          do i=1,size(headers)
2833
2834             !
2835             ! Find entry in METGRID.TBL for this field, if one exists; if it does, then the
2836             !    entry may tell us how to get values for the current field at the missing level
2837             !
2838             do ii=1,num_entries
2839                if (fieldname(ii) == headers(i)%header%field) exit
2840             end do
2841             if (ii <= num_entries) then
2842                call dup(headers(i),new_field)
2843                nullify(new_field%valid_mask)
2844                nullify(new_field%modified_mask)
2845                call fill_field(new_field, ii, output_flags, r_union_levels)
2846             end if
2847
2848          end do
2849
2850          deallocate(union_levels)
2851          deallocate(r_union_levels)
2852          deallocate(headers)
2853
2854          call storage_get_td_headers(headers)
2855
2856          !
2857          ! Now we may need to vertically interpolate to missing values in 3-d fields
2858          !
2859          do i=1,size(headers)
2860
2861             call storage_get_levels(headers(i), field_levels)
2862
2863             ! If this isn't a 3-d array, nothing to do
2864             if (size(field_levels) > 1) then
2865
2866                do k=1,size(field_levels)
2867                   call dup(headers(i),search_field)
2868                   search_field%header%vertical_level = field_levels(k)
2869                   call storage_get_field(search_field,istatus)
2870                   if (istatus == 0) then
2871                      JLOOP: do jx=search_field%header%dim2(1),search_field%header%dim2(2)
2872                         ILOOP: do ix=search_field%header%dim1(1),search_field%header%dim1(2)
2873                            if (.not. bitarray_test(search_field%valid_mask, &
2874                                                    ix-search_field%header%dim1(1)+1, &
2875                                                    jx-search_field%header%dim2(1)+1)) then
2876
2877                               call dup(search_field, lower_field)
2878                               do lower=k-1,1,-1
2879                                  lower_field%header%vertical_level = field_levels(lower)
2880                                  call storage_get_field(lower_field,istatus)
2881                                  if (bitarray_test(lower_field%valid_mask, &
2882                                                    ix-search_field%header%dim1(1)+1, &
2883                                                    jx-search_field%header%dim2(1)+1)) &
2884                                      exit
2885
2886                               end do
2887
2888                               call dup(search_field, upper_field)
2889                               do upper=k+1,size(field_levels)
2890                                  upper_field%header%vertical_level = field_levels(upper)
2891                                  call storage_get_field(upper_field,istatus)
2892                                  if (bitarray_test(upper_field%valid_mask, &
2893                                                    ix-search_field%header%dim1(1)+1, &
2894                                                    jx-search_field%header%dim2(1)+1)) &
2895                                      exit
2896
2897                               end do
2898                               if (upper <= size(field_levels) .and. lower >= 1) then
2899                                  search_field%r_arr(ix,jx) = real(abs(field_levels(upper)-field_levels(k))) &
2900                                                            / real(abs(field_levels(upper)-field_levels(lower))) &
2901                                                            * lower_field%r_arr(ix,jx) &
2902                                                            + real(abs(field_levels(k)-field_levels(lower))) &
2903                                                            / real(abs(field_levels(upper)-field_levels(lower))) &
2904                                                            * upper_field%r_arr(ix,jx)
2905                                  call bitarray_set(search_field%valid_mask, &
2906                                                    ix-search_field%header%dim1(1)+1, &
2907                                                    jx-search_field%header%dim2(1)+1)
2908                               end if
2909                            end if
2910                         end do ILOOP
2911                      end do JLOOP
2912                   else
2913                      call mprintf(.true.,ERROR, &
2914                                   'This is bad, could not get %s at level %i.', &
2915                                   s1=trim(search_field%header%field), i1=field_levels(k))
2916                   end if
2917                end do
2918
2919             end if
2920
2921             deallocate(field_levels)
2922
2923          end do
2924
2925       end if
2926
2927       call list_destroy(temp_levels)
2928       deallocate(all_headers)
2929       deallocate(headers)
2930
2931    end subroutine fill_missing_levels
2932
2933
2934    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2935    ! Name: create_derived_fields
2936    !
2937    ! Purpose:
2938    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2939    subroutine create_derived_fields(arg_gridtype, hdate, xfcst, &
2940                                  we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
2941                                  we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e, &
2942                                  created_this_field, output_flags)
2943
2944       use interp_option_module
2945       use list_module
2946       use module_mergesort
2947       use storage_module
2948
2949       implicit none
2950
2951       ! Arguments
2952       integer, intent(in) :: we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
2953                              we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e
2954       real, intent(in) :: xfcst
2955       logical, dimension(:), intent(inout) :: created_this_field
2956       character (len=1), intent(in) :: arg_gridtype
2957       character (len=24), intent(in) :: hdate
2958       character (len=128), dimension(:), intent(inout) :: output_flags
2959
2960       ! Local variables
2961       integer :: idx, i, j, istatus
2962       type (fg_input) :: field
2963
2964       ! Initialize fg_input structure to store the field
2965       field%header%version = 1
2966       field%header%date = hdate//'        '
2967       field%header%time_dependent = .true.
2968       field%header%mask_field = .false.
2969       field%header%constant_field = .false.
2970       field%header%forecast_hour = xfcst
2971       field%header%fg_source = 'Derived from FG'
2972       field%header%field = ' '
2973       field%header%units = ' '
2974       field%header%description = ' '
2975       field%header%vertical_level = 0
2976       field%header%sr_x = 1
2977       field%header%sr_y = 1
2978       field%header%array_order = 'XY '
2979       field%header%is_wind_grid_rel = .true.
2980       field%header%array_has_missing_values = .false.
2981       nullify(field%r_arr)
2982       nullify(field%valid_mask)
2983       nullify(field%modified_mask)
2984
2985       !
2986       ! Check each entry in METGRID.TBL to see whether it is a derive field
2987       !
2988       do idx=1,num_entries
2989          if (is_derived_field(idx)) then
2990
2991             created_this_field(idx) = .true.
2992
2993             call mprintf(.true.,INFORM,'Going to create the field %s',s1=fieldname(idx))
2994
2995             ! Intialize more fields in storage structure
2996             field%header%field = fieldname(idx)
2997             call get_z_dim_name(fieldname(idx),field%header%vertical_coord)
2998             field%map%stagger = output_stagger(idx)
2999             if (arg_gridtype == 'E') then
3000                field%header%dim1(1) = we_mem_s
3001                field%header%dim1(2) = we_mem_e
3002                field%header%dim2(1) = sn_mem_s
3003                field%header%dim2(2) = sn_mem_e
3004             else if (arg_gridtype == 'C') then
3005                if (output_stagger(idx) == M) then
3006                   field%header%dim1(1) = we_mem_s
3007                   field%header%dim1(2) = we_mem_e
3008                   field%header%dim2(1) = sn_mem_s
3009                   field%header%dim2(2) = sn_mem_e
3010                else if (output_stagger(idx) == U) then
3011                   field%header%dim1(1) = we_mem_stag_s
3012                   field%header%dim1(2) = we_mem_stag_e
3013                   field%header%dim2(1) = sn_mem_s
3014                   field%header%dim2(2) = sn_mem_e
3015                else if (output_stagger(idx) == V) then
3016                   field%header%dim1(1) = we_mem_s
3017                   field%header%dim1(2) = we_mem_e
3018                   field%header%dim2(1) = sn_mem_stag_s
3019                   field%header%dim2(2) = sn_mem_stag_e
3020                end if
3021             end if
3022
3023             call fill_field(field, idx, output_flags)
3024
3025          end if
3026       end do
3027
3028
3029    end subroutine create_derived_fields
3030
3031
3032    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3033    ! Name: fill_field
3034    !
3035    ! Purpose:
3036    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3037    subroutine fill_field(field, idx, output_flags, all_level_list)
3038
3039       use interp_option_module
3040       use list_module
3041       use module_mergesort
3042       use storage_module
3043
3044       implicit none
3045
3046       ! Arguments
3047       integer, intent(in) :: idx
3048       type (fg_input), intent(inout) :: field
3049       character (len=128), dimension(:), intent(inout) :: output_flags
3050       real, dimension(:), intent(in), optional :: all_level_list
3051
3052       ! Local variables
3053       integer :: i, j, istatus, isrclevel
3054       integer, pointer, dimension(:) :: all_list
3055       real :: rfillconst, rlevel, rsrclevel
3056       type (fg_input) :: query_field
3057       type (list_item), pointer, dimension(:) :: keys
3058       character (len=128) :: asrcname
3059       logical :: filled_all_lev
3060
3061       !CWH Initialize local pointer variables
3062       nullify(all_list)
3063       nullify(keys)
3064
3065       filled_all_lev = .false.
3066
3067       !
3068       ! Get a list of all levels to be filled for this field
3069       !
3070       keys => list_get_keys(fill_lev_list(idx))
3071
3072       do i=1,list_length(fill_lev_list(idx))
3073
3074          !
3075          ! First handle a specification for levels "all"
3076          !
3077          if (trim(keys(i)%ckey) == 'all') then
3078
3079             ! We only want to fill all levels if we haven't already filled "all" of them
3080             if (.not. filled_all_lev) then
3081
3082                filled_all_lev = .true.
3083
3084                query_field%header%time_dependent = .true.
3085                query_field%header%field = ' '
3086                nullify(query_field%r_arr)
3087                nullify(query_field%valid_mask)
3088                nullify(query_field%modified_mask)
3089
3090                ! See if we are filling this level with a constant
3091                call get_constant_fill_lev(keys(i)%cvalue, rfillconst, istatus)
3092                if (istatus == 0) then
3093                   if (present(all_level_list)) then
3094                      do j=1,size(all_level_list)
3095                         call create_level(field, real(all_level_list(j)), idx, output_flags, rfillconst=rfillconst)
3096                      end do
3097                   else
3098                      query_field%header%field = level_template(idx)
3099                      nullify(all_list)
3100                      call storage_get_levels(query_field, all_list)
3101                      if (associated(all_list)) then
3102                         do j=1,size(all_list)
3103                            call create_level(field, real(all_list(j)), idx, output_flags, rfillconst=rfillconst)
3104                         end do
3105                         deallocate(all_list)
3106                      end if
3107                   end if
3108
3109                ! Else see if we are filling this level with a constant equal
3110                !   to the value of the level
3111                else if (trim(keys(i)%cvalue) == 'vertical_index') then
3112                   if (present(all_level_list)) then
3113                      do j=1,size(all_level_list)
3114                         call create_level(field, real(all_level_list(j)), idx, output_flags, &
3115                                           rfillconst=real(all_level_list(j)))
3116                      end do
3117                   else
3118                      query_field%header%field = level_template(idx)
3119                      nullify(all_list)
3120                      call storage_get_levels(query_field, all_list)
3121                      if (associated(all_list)) then
3122                         do j=1,size(all_list)
3123                            call create_level(field, real(all_list(j)), idx, output_flags, rfillconst=real(all_list(j)))
3124                         end do
3125                         deallocate(all_list)
3126                      end if
3127                   end if
3128
3129                ! Else, we assume that it is a field from which we are copying levels
3130                else
3131                   if (present(all_level_list)) then
3132                      do j=1,size(all_level_list)
3133                         call create_level(field, real(all_level_list(j)), idx, output_flags, &
3134                                           asrcname=keys(i)%cvalue, rsrclevel=real(all_level_list(j)))
3135                      end do
3136                   else
3137                      query_field%header%field = keys(i)%cvalue  ! Use same levels as source field, not level_template
3138                      nullify(all_list)
3139                      call storage_get_levels(query_field, all_list)
3140                      if (associated(all_list)) then
3141                         do j=1,size(all_list)
3142                            call create_level(field, real(all_list(j)), idx, output_flags, &
3143                                              asrcname=keys(i)%cvalue, rsrclevel=real(all_list(j)))
3144                         end do
3145                         deallocate(all_list)
3146
3147                      else
3148
3149                         ! If the field doesn't have any levels (or does not exist) then we have not
3150                         !   really filled all levels at this point.
3151                         filled_all_lev = .false.
3152                      end if
3153                   end if
3154
3155                end if
3156             end if
3157
3158          !
3159          ! Handle individually specified levels
3160          !
3161          else
3162
3163             read(keys(i)%ckey,*) rlevel
3164
3165             ! See if we are filling this level with a constant
3166             call get_constant_fill_lev(keys(i)%cvalue, rfillconst, istatus)
3167             if (istatus == 0) then
3168                call create_level(field, rlevel, idx, output_flags, rfillconst=rfillconst)
3169
3170             ! Otherwise, we are filling from another level
3171             else
3172                call get_fill_src_level(keys(i)%cvalue, asrcname, isrclevel)
3173                rsrclevel = real(isrclevel)
3174                call create_level(field, rlevel, idx, output_flags, &
3175                                  asrcname=asrcname, rsrclevel=rsrclevel)
3176
3177             end if
3178          end if
3179       end do
3180
3181       if (associated(keys)) deallocate(keys)
3182
3183    end subroutine fill_field
3184
3185
3186    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3187    ! Name: create_level
3188    !
3189    ! Purpose:
3190    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3191    subroutine create_level(field_template, rlevel, idx, output_flags, &
3192                            rfillconst, asrcname, rsrclevel)
3193
3194       use storage_module
3195       use interp_option_module
3196
3197       implicit none
3198
3199       ! Arguments
3200       type (fg_input), intent(inout) :: field_template
3201       real, intent(in) :: rlevel
3202       integer, intent(in) :: idx
3203       character (len=128), dimension(:), intent(inout) :: output_flags
3204       real, intent(in), optional :: rfillconst, rsrclevel
3205       character (len=128), intent(in), optional :: asrcname
3206
3207       ! Local variables
3208       integer :: i, j, istatus
3209       integer :: sm1,em1,sm2,em2
3210       type (fg_input) :: query_field
3211
3212       !
3213       ! Check to make sure optional arguments are sane
3214       !
3215       if (present(rfillconst) .and. (present(asrcname) .or. present(rsrclevel))) then
3216          call mprintf(.true.,ERROR,'A call to create_level() cannot be given specifications '// &
3217                       'for both a constant fill value and a source level.')
3218
3219       else if ((present(asrcname) .and. .not. present(rsrclevel)) .or. &
3220                (.not. present(asrcname) .and. present(rsrclevel))) then
3221          call mprintf(.true.,ERROR,'Neither or both of optional arguments asrcname and '// &
3222                       'rsrclevel must be specified to subroutine create_level().')
3223
3224       else if (.not. present(rfillconst) .and. &
3225                .not. present(asrcname)   .and. &
3226                .not. present(rsrclevel)) then
3227          call mprintf(.true.,ERROR,'A call to create_level() must be given either a specification '// &
3228                       'for a constant fill value or a source level.')
3229       end if
3230
3231       query_field%header%time_dependent = .true.
3232       query_field%header%field = field_template%header%field
3233       query_field%header%vertical_level = rlevel
3234       nullify(query_field%r_arr)
3235       nullify(query_field%valid_mask)
3236       nullify(query_field%modified_mask)
3237
3238       call storage_query_field(query_field, istatus)
3239       if (istatus == 0) then
3240          call mprintf(.true.,INFORM,'%s at level %f already exists; leaving it alone.', &
3241                       s1=field_template%header%field, f1=rlevel)
3242          return
3243       end if
3244
3245       sm1 = field_template%header%dim1(1)
3246       em1 = field_template%header%dim1(2)
3247       sm2 = field_template%header%dim2(1)
3248       em2 = field_template%header%dim2(2)
3249
3250       !
3251       ! Handle constant fill value case
3252       !
3253       if (present(rfillconst)) then
3254
3255          field_template%header%vertical_level = rlevel
3256          allocate(field_template%r_arr(sm1:em1,sm2:em2))
3257          allocate(field_template%valid_mask)
3258          allocate(field_template%modified_mask)
3259          call bitarray_create(field_template%valid_mask, em1-sm1+1, em2-sm2+1)
3260          call bitarray_create(field_template%modified_mask, em1-sm1+1, em2-sm2+1)
3261
3262          field_template%r_arr = rfillconst
3263
3264          do j=sm2,em2
3265             do i=sm1,em1
3266                call bitarray_set(field_template%valid_mask, i-sm1+1, j-sm2+1)
3267             end do
3268          end do
3269
3270          call storage_put_field(field_template)
3271
3272          if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
3273             output_flags(idx) = flag_in_output(idx)
3274          end if
3275
3276       !
3277       ! Handle source field and source level case
3278       !
3279       else if (present(asrcname) .and. present(rsrclevel)) then
3280
3281          query_field%header%field = ' '
3282          query_field%header%field = asrcname
3283          query_field%header%vertical_level = rsrclevel
3284
3285          ! Check to see whether the requested source field exists at the requested level
3286          call storage_query_field(query_field, istatus)
3287
3288          if (istatus == 0) then
3289
3290             ! Read in requested field at requested level
3291             call storage_get_field(query_field, istatus)
3292             if ((query_field%header%dim1(1) /= field_template%header%dim1(1)) .or. &
3293                 (query_field%header%dim1(2) /= field_template%header%dim1(2)) .or. &
3294                 (query_field%header%dim2(1) /= field_template%header%dim2(1)) .or. &
3295                 (query_field%header%dim2(2) /= field_template%header%dim2(2))) then
3296                call mprintf(.true.,ERROR,'Dimensions for %s do not match those of %s. This is '// &
3297                             'probably because the staggerings of the fields do not match.', &
3298                             s1=query_field%header%field, s2=field_template%header%field)
3299             end if
3300
3301             field_template%header%vertical_level = rlevel
3302             allocate(field_template%r_arr(sm1:em1,sm2:em2))
3303             allocate(field_template%valid_mask)
3304             allocate(field_template%modified_mask)
3305             call bitarray_create(field_template%valid_mask, em1-sm1+1, em2-sm2+1)
3306             call bitarray_create(field_template%modified_mask, em1-sm1+1, em2-sm2+1)
3307
3308             field_template%r_arr = query_field%r_arr
3309
3310             ! We should retain information about which points in the field are valid
3311             do j=sm2,em2
3312                do i=sm1,em1
3313                   if (bitarray_test(query_field%valid_mask, i-sm1+1, j-sm2+1)) then
3314                      call bitarray_set(field_template%valid_mask, i-sm1+1, j-sm2+1)
3315                   end if
3316                end do
3317             end do
3318
3319             call storage_put_field(field_template)
3320
3321             if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
3322                output_flags(idx) = flag_in_output(idx)
3323             end if
3324
3325          else
3326             call mprintf(.true.,INFORM,'Couldn''t find %s at level %f to fill level %f of %s.', &
3327                          s1=asrcname,f1=rsrclevel,f2=rlevel,s2=field_template%header%field)
3328          end if
3329
3330       end if
3331
3332    end subroutine create_level
3333
3334
3335    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3336    ! Name: accum_continuous
3337    !
3338    ! Purpose: Sum up all of the source data points whose nearest neighbor in the
3339    !   model grid is the specified model grid point.
3340    !
3341    ! NOTE: When processing the source tile, those source points that are
3342    !   closest to a different model grid point will be added to the totals for
3343    !   such grid points; thus, an entire source tile will be processed at a time.
3344    !   This routine really processes for all model grid points that are
3345    !   within a source tile, and not just for a single grid point.
3346    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3347    subroutine accum_continuous(src_array, &
3348                                src_min_x, src_max_x, src_min_y, src_max_y, src_min_z, src_max_z, bdr_width, &
3349                                dst_array, n, &
3350                                start_i, end_i, start_j, end_j, start_k, end_k, &
3351                                istagger, &
3352                                new_pts, msgval, maskval, mask_relational, mask_array, sr_x, sr_y)
3353
3354       use bitarray_module
3355       use misc_definitions_module
3356
3357       implicit none
3358
3359       ! Arguments
3360       integer, intent(in) :: start_i, end_i, start_j, end_j, start_k, end_k, istagger, &
3361                              src_min_x, src_max_x, src_min_y, src_max_y, src_min_z, src_max_z, bdr_width
3362       real, intent(in) :: maskval, msgval
3363       real, dimension(src_min_x:src_max_x, src_min_y:src_max_y, src_min_z:src_max_z), intent(in) :: src_array
3364       real, dimension(start_i:end_i, start_j:end_j, start_k:end_k), intent(inout) :: dst_array, n
3365       real, dimension(src_min_x:src_max_x, src_min_y:src_max_y), intent(in), optional :: mask_array
3366       integer, intent(in), optional :: sr_x, sr_y
3367       type (bitarray), intent(inout) :: new_pts
3368       character(len=1), intent(in), optional :: mask_relational
3369
3370       ! Local variables
3371       integer :: i, j
3372       integer, pointer, dimension(:,:,:) :: where_maps_to
3373       real :: rsr_x, rsr_y
3374
3375       rsr_x = 1.0
3376       rsr_y = 1.0
3377       if (present(sr_x)) rsr_x = real(sr_x)
3378       if (present(sr_y)) rsr_y = real(sr_y)
3379
3380       allocate(where_maps_to(src_min_x:src_max_x,src_min_y:src_max_y,2))
3381       do i=src_min_x,src_max_x
3382          do j=src_min_y,src_max_y
3383             where_maps_to(i,j,1) = NOT_PROCESSED
3384          end do
3385       end do
3386
3387       call process_continuous_block(src_array, where_maps_to, &
3388                                src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
3389                                src_min_x+bdr_width, src_min_y, src_min_z, &
3390                                src_max_x-bdr_width, src_max_y, src_max_z, &
3391                                dst_array, n, start_i, end_i, start_j, end_j, start_k, end_k, &
3392                                istagger, &
3393                                new_pts, rsr_x, rsr_y, msgval, maskval, mask_relational, mask_array)
3394
3395       deallocate(where_maps_to)
3396
3397    end subroutine accum_continuous
3398
3399
3400    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3401    ! Name: process_continuous_block
3402    !
3403    ! Purpose: To recursively process a subarray of continuous data, adding the
3404    !   points in a block to the sum for their nearest grid point. The nearest
3405    !   neighbor may be estimated in some cases; for example, if the four corners
3406    !   of a subarray all have the same nearest grid point, all elements in the
3407    !   subarray are added to that grid point.
3408    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3409    recursive subroutine process_continuous_block(tile_array, where_maps_to, &
3410                                    src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
3411                                    min_i, min_j, min_k, max_i, max_j, max_k, &
3412                                    dst_array, n, &
3413                                    start_x, end_x, start_y, end_y, start_z, end_z, &
3414                                    istagger, &
3415                                    new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
3416
3417       use bitarray_module
3418       use llxy_module
3419       use misc_definitions_module
3420
3421       implicit none
3422
3423       ! Arguments
3424       integer, intent(in) :: min_i, min_j, min_k, max_i, max_j, max_k, &
3425                              src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
3426                              start_x, end_x, start_y, end_y, start_z, end_z, istagger
3427       integer, dimension(src_min_x:src_max_x,src_min_y:src_max_y,2), intent(inout) :: where_maps_to
3428       real, intent(in) :: sr_x, sr_y, maskval, msgval
3429       real, dimension(src_min_x:src_max_x,src_min_y:src_max_y,src_min_z:src_max_z), intent(in) :: tile_array
3430       real, dimension(src_min_x:src_max_x,src_min_y:src_max_y), intent(in), optional :: mask_array
3431       real, dimension(start_x:end_x,start_y:end_y,start_z:end_z), intent(inout) :: dst_array, n
3432       type (bitarray), intent(inout) :: new_pts
3433       character(len=1), intent(in), optional :: mask_relational
3434
3435       ! Local variables
3436       integer :: orig_selected_domain, x_dest, y_dest, i, j, k, center_i, center_j
3437       real :: lat_corner, lon_corner, rx, ry
3438
3439       ! Compute the model grid point that the corners of the rectangle to be
3440       !   processed map to
3441       ! Lower-left corner
3442       if (where_maps_to(min_i,min_j,1) == NOT_PROCESSED) then
3443          orig_selected_domain = iget_selected_domain()
3444          call select_domain(SOURCE_PROJ)
3445          call xytoll(real(min_i), real(min_j), lat_corner, lon_corner, istagger)
3446          call select_domain(1)
3447          call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
3448          rx = (rx - 0.5)*sr_x + 0.5
3449          ry = (ry - 0.5)*sr_y + 0.5
3450          call select_domain(orig_selected_domain)
3451          if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
3452              real(start_y) <= ry .and. ry <= real(end_y)) then
3453             where_maps_to(min_i,min_j,1) = nint(rx)
3454             where_maps_to(min_i,min_j,2) = nint(ry)
3455          else
3456             where_maps_to(min_i,min_j,1) = OUTSIDE_DOMAIN
3457          end if
3458       end if
3459
3460       ! Upper-left corner
3461       if (where_maps_to(min_i,max_j,1) == NOT_PROCESSED) then
3462          orig_selected_domain = iget_selected_domain()
3463          call select_domain(SOURCE_PROJ)
3464          call xytoll(real(min_i), real(max_j), lat_corner, lon_corner, istagger)
3465          call select_domain(1)
3466          call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
3467          rx = (rx - 0.5)*sr_x + 0.5
3468          ry = (ry - 0.5)*sr_y + 0.5
3469          call select_domain(orig_selected_domain)
3470          if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
3471              real(start_y) <= ry .and. ry <= real(end_y)) then
3472             where_maps_to(min_i,max_j,1) = nint(rx)
3473             where_maps_to(min_i,max_j,2) = nint(ry)
3474          else
3475             where_maps_to(min_i,max_j,1) = OUTSIDE_DOMAIN
3476          end if
3477       end if
3478
3479       ! Upper-right corner
3480       if (where_maps_to(max_i,max_j,1) == NOT_PROCESSED) then
3481          orig_selected_domain = iget_selected_domain()
3482          call select_domain(SOURCE_PROJ)
3483          call xytoll(real(max_i), real(max_j), lat_corner, lon_corner, istagger)
3484          call select_domain(1)
3485          call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
3486          rx = (rx - 0.5)*sr_x + 0.5
3487          ry = (ry - 0.5)*sr_y + 0.5
3488          call select_domain(orig_selected_domain)
3489          if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
3490              real(start_y) <= ry .and. ry <= real(end_y)) then
3491             where_maps_to(max_i,max_j,1) = nint(rx)
3492             where_maps_to(max_i,max_j,2) = nint(ry)
3493          else
3494             where_maps_to(max_i,max_j,1) = OUTSIDE_DOMAIN
3495          end if
3496       end if
3497
3498       ! Lower-right corner
3499       if (where_maps_to(max_i,min_j,1) == NOT_PROCESSED) then
3500          orig_selected_domain = iget_selected_domain()
3501          call select_domain(SOURCE_PROJ)
3502          call xytoll(real(max_i), real(min_j), lat_corner, lon_corner, istagger)
3503          call select_domain(1)
3504          call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
3505          rx = (rx - 0.5)*sr_x + 0.5
3506          ry = (ry - 0.5)*sr_y + 0.5
3507          call select_domain(orig_selected_domain)
3508          if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
3509              real(start_y) <= ry .and. ry <= real(end_y)) then
3510             where_maps_to(max_i,min_j,1) = nint(rx)
3511             where_maps_to(max_i,min_j,2) = nint(ry)
3512          else
3513             where_maps_to(max_i,min_j,1) = OUTSIDE_DOMAIN
3514          end if
3515       end if
3516
3517       ! If all four corners map to same model grid point, accumulate the
3518       !   entire rectangle
3519       if (where_maps_to(min_i,min_j,1) == where_maps_to(min_i,max_j,1) .and. &
3520           where_maps_to(min_i,min_j,1) == where_maps_to(max_i,max_j,1) .and. &
3521           where_maps_to(min_i,min_j,1) == where_maps_to(max_i,min_j,1) .and. &
3522           where_maps_to(min_i,min_j,2) == where_maps_to(min_i,max_j,2) .and. &
3523           where_maps_to(min_i,min_j,2) == where_maps_to(max_i,max_j,2) .and. &
3524           where_maps_to(min_i,min_j,2) == where_maps_to(max_i,min_j,2) .and. &
3525           where_maps_to(min_i,min_j,1) /= OUTSIDE_DOMAIN) then
3526          x_dest = where_maps_to(min_i,min_j,1)
3527          y_dest = where_maps_to(min_i,min_j,2)
3528
3529          ! If this grid point was already given a value from higher-priority source data,
3530          !   there is nothing to do.
3531 !         if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
3532
3533             ! If this grid point has never been given a value by this level of source data,
3534             !   initialize the point
3535             if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
3536                do k=min_k,max_k
3537                   dst_array(x_dest,y_dest,k) = 0.
3538                end do
3539             end if
3540
3541             ! Sum all the points whose nearest neighbor is this grid point
3542             if (present(mask_array) .and. present(mask_relational)) then
3543                do i=min_i,max_i
3544                   do j=min_j,max_j
3545                      do k=min_k,max_k
3546                         ! Ignore masked/missing values in the source data
3547                         if (tile_array(i,j,k) /= msgval) then
3548                            if (mask_relational == ' ' .and. mask_array(i,j) /= maskval) then
3549                               dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3550                               n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3551                               call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3552                            else if (mask_relational == '<' .and. mask_array(i,j) >= maskval) then
3553                               dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3554                               n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3555                               call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3556                            else if (mask_relational == '>' .and. mask_array(i,j) <= maskval) then
3557                               dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3558                               n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3559                               call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3560                            end if
3561                         end if
3562                      end do
3563                   end do
3564                end do
3565             else if (present(mask_array)) then
3566                do i=min_i,max_i
3567                   do j=min_j,max_j
3568                      do k=min_k,max_k
3569                         ! Ignore masked/missing values in the source data
3570                         if ((tile_array(i,j,k) /= msgval) .and. &
3571                             (mask_array(i,j) /= maskval)) then
3572                            dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3573                            n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3574                            call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3575                         end if
3576                      end do
3577                   end do
3578                end do
3579             else
3580                do i=min_i,max_i
3581                   do j=min_j,max_j
3582                      do k=min_k,max_k
3583                         ! Ignore masked/missing values in the source data
3584                         if ((tile_array(i,j,k) /= msgval)) then
3585                            dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3586                            n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3587                            call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3588                         end if
3589                      end do
3590                   end do
3591                end do
3592             end if
3593
3594 !         end if
3595
3596       ! Rectangle is a square of four points, and we can simply deal with each of the points
3597       else if (((max_i - min_i + 1) <= 2) .and. ((max_j - min_j + 1) <= 2)) then
3598          do i=min_i,max_i
3599             do j=min_j,max_j
3600                x_dest = where_maps_to(i,j,1)
3601                y_dest = where_maps_to(i,j,2)
3602
3603                if (x_dest /= OUTSIDE_DOMAIN) then
3604
3605 !                  if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
3606                      if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
3607                         do k=min_k,max_k
3608                            dst_array(x_dest,y_dest,k) = 0.
3609                         end do
3610                      end if
3611
3612                      if (present(mask_array) .and. present(mask_relational)) then
3613                         do k=min_k,max_k
3614                            ! Ignore masked/missing values
3615                            if (tile_array(i,j,k) /= msgval) then
3616                               if (mask_relational == ' ' .and. mask_array(i,j) /= maskval) then
3617                                  dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3618                                  n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3619                                  call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3620                               else if (mask_relational == '<' .and. mask_array(i,j) >= maskval) then
3621                                  dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3622                                  n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3623                                  call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3624                               else if (mask_relational == '>' .and. mask_array(i,j) <= maskval) then
3625                                  dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3626                                  n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3627                                  call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3628                               end if
3629                            end if
3630                         end do
3631                      else if (present(mask_array)) then
3632                         do k=min_k,max_k
3633                            ! Ignore masked/missing values
3634                            if ((tile_array(i,j,k) /= msgval) .and. &
3635                                 (mask_array(i,j) /= maskval)) then
3636                               dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3637                               n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3638                               call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3639                            end if
3640                         end do
3641                      else
3642                         do k=min_k,max_k
3643                            ! Ignore masked/missing values
3644                            if ((tile_array(i,j,k) /= msgval)) then
3645                               dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3646                               n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3647                               call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3648                            end if
3649                         end do
3650                      end if
3651 !                  end if
3652
3653                end if
3654             end do
3655          end do
3656
3657       ! Not all corners map to the same grid point, and the rectangle contains more than
3658       !   four points
3659       else
3660          center_i = (max_i + min_i)/2
3661          center_j = (max_j + min_j)/2
3662
3663          ! Recursively process lower-left rectangle
3664          call process_continuous_block(tile_array, where_maps_to, &
3665                     src_min_x, src_min_y, src_min_z, &
3666                     src_max_x, src_max_y, src_max_z, &
3667                     min_i, min_j, min_k, &
3668                     center_i, center_j, max_k, &
3669                     dst_array, n, &
3670                     start_x, end_x, start_y, end_y, start_z, end_z, &
3671                     istagger, &
3672                     new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
3673
3674          if (center_i < max_i) then
3675             ! Recursively process lower-right rectangle
3676             call process_continuous_block(tile_array, where_maps_to, &
3677                        src_min_x, src_min_y, src_min_z, &
3678                        src_max_x, src_max_y, src_max_z, &
3679                        center_i+1, min_j, min_k, max_i, &
3680                        center_j, max_k, &
3681                        dst_array, n, &
3682                        start_x, end_x, start_y, &
3683                        end_y, start_z, end_z, &
3684                        istagger, &
3685                        new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
3686          end if
3687
3688          if (center_j < max_j) then
3689             ! Recursively process upper-left rectangle
3690             call process_continuous_block(tile_array, where_maps_to, &
3691                        src_min_x, src_min_y, src_min_z, &
3692                        src_max_x, src_max_y, src_max_z, &
3693                        min_i, center_j+1, min_k, center_i, &
3694                        max_j, max_k, &
3695                        dst_array, n, &
3696                        start_x, end_x, start_y, &
3697                        end_y, start_z, end_z, &
3698                        istagger, &
3699                        new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
3700          end if
3701
3702          if (center_i < max_i .and. center_j < max_j) then
3703             ! Recursively process upper-right rectangle
3704             call process_continuous_block(tile_array, where_maps_to, &
3705                        src_min_x, src_min_y, src_min_z, &
3706                        src_max_x, src_max_y, src_max_z, &
3707                        center_i+1, center_j+1, min_k, max_i, &
3708                        max_j, max_k, &
3709                        dst_array, n, &
3710                        start_x, end_x, start_y, &
3711                        end_y, start_z, end_z, &
3712                        istagger, &
3713                        new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
3714          end if
3715       end if
3716
3717    end subroutine process_continuous_block
3718
3719 end module process_domain_module