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) :: mpas_source_mesh
   8    type (target_mesh_type) :: wrf_target_mesh_m, wrf_target_mesh_u, wrf_target_mesh_v
   9    type (remap_info_type) :: 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       if (.not. do_const_processing) then
 855          call derive_mpas_fields()
 856       end if
 857
 858       !
 859       ! Check that every mandatory field was found in input data
 860       !
 861       do i=1,num_entries
 862          if (is_mandatory(i) .and. .not. got_this_field(i)) then
 863             call mprintf(.true.,ERROR,'The mandatory field %s was not found in any input data.',s1=fieldname(i))
 864          end if
 865       end do
 866
 867       !
 868       ! Before we begin to write fields, if debug_level is set high enough, we
 869       !    write a table of which fields are available at which levels to the
 870       !    metgrid.log file, and then we check to see if any fields are not
 871       !    completely covered with data.
 872       !
 873       call storage_print_fields()
 874       call find_missing_values()
 875
 876       !
 877       ! All of the processing is now done for this time period for this domain;
 878       !   now we simply output every field from the storage module.
 879       !
 880
 881       title = 'OUTPUT FROM METGRID V3.9'
 882
 883       ! Initialize the output module for this domain and time
 884       call mprintf(.true.,LOGFILE,'Initializing output module.')
 885       output_date = temp_date
 886       if ( .not. nocolons ) then
 887          if (len_trim(temp_date) == 13) then
 888             output_date(14:19) = ':00:00'
 889          else if (len_trim(temp_date) == 16) then
 890             output_date(17:19) = ':00'
 891          end if
 892       else
 893          if (len_trim(temp_date) == 13) then
 894             output_date(14:19) = '_00_00'
 895          else if (len_trim(temp_date) == 16) then
 896             output_date(17:19) = '_00'
 897          end if
 898       endif
 899
 900       call output_init(n, title, output_date, gridtype, dyn_opt, &
 901                        corner_lats, corner_lons, &
 902                        we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
 903                        we_patch_s,  we_patch_e,  sn_patch_s,  sn_patch_e, &
 904                        we_mem_s,    we_mem_e,    sn_mem_s,    sn_mem_e, &
 905                        extra_col, extra_row)
 906
 907       call get_bottom_top_dim(bottom_top_dim)
 908
 909       ! Add in a flag to tell real that we have seven new msf fields
 910       nmet_flags = num_entries + 1
 911       if (associated(geogrid_flags)) nmet_flags = nmet_flags + size(geogrid_flags)
 912       allocate(met_flags(nmet_flags))
 913       do i=1,num_entries
 914          met_flags(i) = output_flags(i)
 915       end do
 916       if (gridtype == 'C') then
 917          met_flags(num_entries+1) = 'FLAG_MF_XY'
 918       else
 919          met_flags(num_entries+1) = ' '
 920       end if
 921       if (associated(geogrid_flags)) then
 922          do i=1,size(geogrid_flags)
 923             met_flags(num_entries+1+i) = geogrid_flags(i)
 924          end do
 925       end if
 926
 927       ! Find out how many soil levels or layers we have; this assumes a field named ST
 928       field % header % field = 'ST'
 929       nullify(soil_levels)
 930       call storage_get_levels(field, soil_levels)
 931
 932       if (.not. associated(soil_levels)) then
 933          field % header % field = 'SOILT'
 934          nullify(soil_levels)
 935          call storage_get_levels(field, soil_levels)
 936       end if
 937
 938       if (.not. associated(soil_levels)) then
 939          field % header % field = 'STC_WPS'
 940          nullify(soil_levels)
 941          call storage_get_levels(field, soil_levels)
 942       end if
 943
 944       if (associated(soil_levels)) then
 945          num_metgrid_soil_levs = size(soil_levels)
 946          deallocate(soil_levels)
 947       else
 948          num_metgrid_soil_levs = 0
 949       end if
 950
 951       ! First write out global attributes
 952       call mprintf(.true.,LOGFILE,'Writing global attributes to output.')
 953       call write_global_attrs(title, output_date, gridtype, dyn_opt, west_east_dim,          &
 954                               south_north_dim, bottom_top_dim,                               &
 955                               we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e,      &
 956                               sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e,      &
 957                               map_proj, mminlu, num_land_cat,                                &
 958                               is_water, is_lake, is_ice, is_urban, i_soilwater,              &
 959                               grid_id, parent_id, i_parent_start,                            &
 960                               j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy,    &
 961                               cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
 962                               pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y,           &
 963                               corner_lats, corner_lons, num_metgrid_soil_levs,               &
 964                               met_flags, nmet_flags, process_bdy_width)
 965
 966       deallocate(met_flags)
 967
 968       call reset_next_field()
 969
 970       istatus = 0
 971
 972       ! Now loop over all output fields, writing each to the output module
 973       do while (istatus == 0)
 974          call get_next_output_field(cname, real_array, &
 975                                     sm1, em1, sm2, em2, sm3, em3, istatus)
 976          if (istatus == 0) then
 977
 978             call mprintf(.true.,LOGFILE,'Writing field %s to output.',s1=cname)
 979             call write_field(sm1, em1, sm2, em2, sm3, em3, &
 980                              cname, output_date, real_array)
 981             deallocate(real_array)
 982
 983          end if
 984
 985       end do
 986
 987       call mprintf(.true.,LOGFILE,'Closing output file.')
 988       call output_close()
 989
 990       ! Free up memory used by met fields for this valid time
 991       call storage_delete_all_td()
 992
 993    end subroutine process_single_met_time
 994
 995
 996    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 997    ! Name: process_intermediate_fields
 998    !
 999    ! Purpose:
1000    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1001    subroutine process_intermediate_fields(input_name, do_const_processing, temp_date, fg_data, got_this_field, &
1002                                           landmask, process_bdy_width, &
1003                                           u_field, v_field, &
1004                                           dom_dx, dom_dy, &
1005                                           xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v, &
1006                                           output_flags, west_east_dim, south_north_dim, &
1007                                           we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1008                                           we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
1009                                           sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
1010                                           we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
1011                                           sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e )
1012
1013       use bitarray_module
1014       use gridinfo_module
1015       use interp_module
1016       use interp_option_module
1017       use llxy_module
1018       use misc_definitions_module
1019       use module_debug
1020       use output_module
1021       use parallel_module
1022       use read_met_module
1023       use rotate_winds_module
1024       use storage_module
1025
1026       implicit none
1027
1028 ! BUG: Move this constant to misc_definitions_module?
1029 integer, parameter :: BDR_WIDTH = 3
1030
1031       character (len=*), intent(inout) :: input_name
1032       logical, intent(in) :: do_const_processing
1033       character (len=*), intent(in) :: temp_date
1034       type (met_data), intent(inout) :: fg_data
1035       logical, dimension(:), intent(inout) :: got_this_field
1036       real, pointer, dimension(:,:) :: landmask
1037       integer, intent(in) :: process_bdy_width
1038       type (fg_input), intent(inout) :: u_field, v_field
1039       character (len=128), dimension(:), intent(inout) :: output_flags
1040       real, intent(in) :: dom_dx, dom_dy
1041       real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
1042       integer, intent(in) :: west_east_dim, south_north_dim, &
1043                  we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1044                  we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
1045                  sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
1046                  we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
1047                  sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e
1048
1049       integer :: istatus
1050       integer :: idx, idxt, u_idx
1051       integer :: iqstatus
1052       real :: threshold
1053       real, pointer, dimension(:,:) :: halo_slab => null()
1054       integer, pointer, dimension(:) :: u_levels => null()
1055       integer, pointer, dimension(:) :: v_levels => null()
1056       integer :: bdr_wdth
1057       logical :: do_gcell_interp
1058       type (fg_input) :: field
1059
1060
1061       ! Do a first pass through this fg source to get all mask fields used
1062       !   during interpolation
1063       call get_interp_masks(trim(input_name), do_const_processing, temp_date)
1064
1065       istatus = 0
1066
1067       ! Initialize the module for reading in the met fields
1068       call read_met_init(trim(input_name), do_const_processing, temp_date, istatus)
1069
1070       if (istatus == 0) then
1071
1072          ! Process all fields and levels from the current file; read_next_met_field()
1073          !   will return a non-zero status when there are no more fields to be read.
1074          do while (istatus == 0)
1075
1076             call read_next_met_field(fg_data, istatus)
1077
1078             if (istatus == 0) then
1079
1080                ! Find index into fieldname, interp_method, masked, and fill_missing
1081                !   of the current field
1082                idxt = num_entries + 1
1083                do idx=1,num_entries
1084                   if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
1085                       (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
1086
1087                      got_this_field(idx) = .true.
1088
1089                      if (index(input_name,trim(from_input(idx))) /= 0 .or. &
1090                         (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1091                         idxt = idx
1092                      end if
1093
1094                   end if
1095                end do
1096                idx = idxt
1097                if (idx > num_entries) idx = num_entries ! The last entry is a default
1098
1099                ! Do we need to rename this field?
1100                if (output_name(idx) /= ' ') then
1101                   fg_data%field = output_name(idx)(1:9)
1102
1103                   idxt = num_entries + 1
1104                   do idx=1,num_entries
1105                      if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
1106                          (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
1107
1108                         got_this_field(idx) = .true.
1109
1110                         if (index(input_name,trim(from_input(idx))) /= 0 .or. &
1111                            (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1112                            idxt = idx
1113                         end if
1114
1115                      end if
1116                   end do
1117                   idx = idxt
1118                   if (idx > num_entries) idx = num_entries ! The last entry is a default
1119                end if
1120
1121                ! Do a simple check to see whether this is a global dataset
1122                ! Note that we do not currently support regional Gaussian grids
1123                if ((fg_data%iproj == PROJ_LATLON .and. abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) &
1124                     .or. (fg_data%iproj == PROJ_GAUSS)) then
1125                   bdr_wdth = BDR_WIDTH
1126                   allocate(halo_slab(1-BDR_WIDTH:fg_data%nx+BDR_WIDTH,1:fg_data%ny))
1127
1128                   halo_slab(1:fg_data%nx,                      1:fg_data%ny) = &
1129                             fg_data%slab(1:fg_data%nx,              1:fg_data%ny)
1130
1131                   halo_slab(1-BDR_WIDTH:0,                     1:fg_data%ny) = &
1132                             fg_data%slab(fg_data%nx-BDR_WIDTH+1:fg_data%nx, 1:fg_data%ny)
1133
1134                   halo_slab(fg_data%nx+1:fg_data%nx+BDR_WIDTH, 1:fg_data%ny) = &
1135                             fg_data%slab(1:BDR_WIDTH,       1:fg_data%ny)
1136
1137                   deallocate(fg_data%slab)
1138                else
1139                   bdr_wdth = 0
1140                   halo_slab => fg_data%slab
1141                   nullify(fg_data%slab)
1142                end if
1143
1144                call mprintf(.true.,LOGFILE,'Processing %s at level %f.',s1=fg_data%field,f1=fg_data%xlvl)
1145
1146                call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, &
1147                                            fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, &
1148                                            fg_data%deltalon, fg_data%starti, fg_data%startj, &
1149                                            fg_data%startlat, fg_data%startlon,  &
1150                                            fg_data%pole_lat, fg_data%pole_lon,        &
1151                                            fg_data%centerlat, fg_data%centerlon,      &
1152                                            real(fg_data%nx+1)/2., real(fg_data%ny+1)/2.,  &
1153                                            earth_radius=fg_data%earth_radius*1000.)
1154
1155                ! Initialize fg_input structure to store the field
1156                field%header%version = 1
1157                field%header%date = fg_data%hdate//'        '
1158                if (do_const_processing) then
1159                   field%header%time_dependent = .false.
1160                   field%header%constant_field = .true.
1161                else
1162                   field%header%time_dependent = .true.
1163                   field%header%constant_field = .false.
1164                end if
1165                field%header%forecast_hour = fg_data%xfcst
1166                field%header%fg_source = 'FG'
1167                field%header%field = ' '
1168                field%header%field(1:9) = fg_data%field
1169                field%header%units = ' '
1170                field%header%units(1:25) = fg_data%units
1171                field%header%description = ' '
1172                field%header%description(1:46) = fg_data%desc
1173                call get_z_dim_name(fg_data%field,field%header%vertical_coord)
1174                field%header%vertical_level = nint(fg_data%xlvl)
1175                field%header%sr_x = 1
1176                field%header%sr_y = 1
1177                field%header%array_order = 'XY '
1178                field%header%is_wind_grid_rel = fg_data%is_wind_grid_rel
1179                field%header%array_has_missing_values = .false.
1180                nullify(field%r_arr)
1181                nullify(field%valid_mask)
1182                nullify(field%modified_mask)
1183
1184                if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
1185                   output_flags(idx) = flag_in_output(idx)
1186                end if
1187
1188                ! If we should not output this field, just list it as a mask field
1189                if (output_this_field(idx)) then
1190                   field%header%mask_field = .false.
1191                else
1192                   field%header%mask_field = .true.
1193                end if
1194
1195                !
1196                ! Before actually doing any interpolation to the model grid, we must check
1197                !    whether we will be using the average_gcell interpolator that averages all
1198                !    source points in each model grid cell
1199                !
1200                do_gcell_interp = .false.
1201                if (index(interp_method(idx),'average_gcell') /= 0) then
1202
1203                   call get_gcell_threshold(interp_method(idx), threshold, istatus)
1204                   if (istatus == 0) then
1205                      if (fg_data%dx == 0. .and. fg_data%dy == 0. .and. &
1206                          fg_data%deltalat /= 0. .and. fg_data%deltalon /= 0.) then
1207                         fg_data%dx = abs(fg_data%deltalon)
1208                         fg_data%dy = abs(fg_data%deltalat)
1209                      else
1210 ! BUG: Need to more correctly handle dx/dy in meters.
1211                         fg_data%dx = fg_data%dx / 111000.  ! Convert meters to approximate degrees
1212                         fg_data%dy = fg_data%dy / 111000.
1213                      end if
1214                      if (gridtype == 'C') then
1215                         if (threshold*max(fg_data%dx,fg_data%dy)*111. <= max(dom_dx,dom_dy)/1000.) &
1216                            do_gcell_interp = .true.
1217                      else if (gridtype == 'E') then
1218                         if (threshold*max(fg_data%dx,fg_data%dy) <= max(dom_dx,dom_dy)) &
1219                            do_gcell_interp = .true.
1220                      end if
1221                   end if
1222                end if
1223
1224                ! Interpolate to U staggering
1225                if (output_stagger(idx) == U) then
1226
1227                   call storage_query_field(field, iqstatus)
1228                   if (iqstatus == 0) then
1229                      call storage_get_field(field, iqstatus)
1230                      call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
1231                                   ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
1232                      if (associated(field%modified_mask)) then
1233                         call bitarray_destroy(field%modified_mask)
1234                         nullify(field%modified_mask)
1235                      end if
1236                   else
1237                      allocate(field%valid_mask)
1238                      call bitarray_create(field%valid_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
1239                   end if
1240
1241                   ! Save a copy of the fg_input structure for the U field so that we can find it later
1242                   if (is_u_field(idx)) call dup(field, u_field)
1243
1244                   allocate(field%modified_mask)
1245                   call bitarray_create(field%modified_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
1246
1247                   if (do_const_processing .or. field%header%time_dependent) then
1248                      call interp_met_field(input_name, fg_data%field, U, M, &
1249                                   field, xlat_u, xlon_u, we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, &
1250                                   we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1251                                   halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
1252                                   field%modified_mask, process_bdy_width)
1253                   else
1254                      call mprintf(.true.,INFORM,' - already processed this field from constant file.')
1255                   end if
1256
1257                ! Interpolate to V staggering
1258                else if (output_stagger(idx) == V) then
1259
1260                   call storage_query_field(field, iqstatus)
1261                   if (iqstatus == 0) then
1262                      call storage_get_field(field, iqstatus)
1263                      call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
1264                                   ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
1265                      if (associated(field%modified_mask)) then
1266                         call bitarray_destroy(field%modified_mask)
1267                         nullify(field%modified_mask)
1268                      end if
1269                   else
1270                      allocate(field%valid_mask)
1271                      call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
1272                   end if
1273
1274                   ! Save a copy of the fg_input structure for the V field so that we can find it later
1275                   if (is_v_field(idx)) call dup(field, v_field)
1276
1277                   allocate(field%modified_mask)
1278                   call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
1279
1280                   if (do_const_processing .or. field%header%time_dependent) then
1281                      call interp_met_field(input_name, fg_data%field, V, M, &
1282                                   field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
1283                                   we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1284                                   halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
1285                                   field%modified_mask, process_bdy_width)
1286                   else
1287                      call mprintf(.true.,INFORM,' - already processed this field from constant file.')
1288                   end if
1289
1290                ! Interpolate to VV staggering
1291                else if (output_stagger(idx) == VV) then
1292
1293                   call storage_query_field(field, iqstatus)
1294                   if (iqstatus == 0) then
1295                      call storage_get_field(field, iqstatus)
1296                      call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
1297                                   ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
1298                      if (associated(field%modified_mask)) then
1299                         call bitarray_destroy(field%modified_mask)
1300                         nullify(field%modified_mask)
1301                      end if
1302                   else
1303                      allocate(field%valid_mask)
1304                      call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
1305                   end if
1306
1307                   ! Save a copy of the fg_input structure for the U field so that we can find it later
1308                   if (is_u_field(idx)) call dup(field, u_field)
1309
1310                   ! Save a copy of the fg_input structure for the V field so that we can find it later
1311                   if (is_v_field(idx)) call dup(field, v_field)
1312
1313                   allocate(field%modified_mask)
1314                   call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
1315
1316                   if (do_const_processing .or. field%header%time_dependent) then
1317                      call interp_met_field(input_name, fg_data%field, VV, M, &
1318                                   field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1319                                   we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1320                                   halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
1321                                   field%modified_mask, process_bdy_width)
1322                   else
1323                      call mprintf(.true.,INFORM,' - already processed this field from constant file.')
1324                   end if
1325
1326                ! All other fields interpolated to M staggering for C grid, H staggering for E grid
1327                else
1328
1329                   call storage_query_field(field, iqstatus)
1330                   if (iqstatus == 0) then
1331                      call storage_get_field(field, iqstatus)
1332                      call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
1333                                   ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
1334                      if (associated(field%modified_mask)) then
1335                         call bitarray_destroy(field%modified_mask)
1336                         nullify(field%modified_mask)
1337                      end if
1338                   else
1339                      allocate(field%valid_mask)
1340                      call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
1341                   end if
1342
1343                   allocate(field%modified_mask)
1344                   call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
1345
1346                   if (do_const_processing .or. field%header%time_dependent) then
1347                      if (gridtype == 'C') then
1348                         call interp_met_field(input_name, fg_data%field, M, M, &
1349                                      field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1350                                      we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1351                                      halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
1352                                      field%modified_mask, process_bdy_width, landmask)
1353
1354                      else if (gridtype == 'E') then
1355                         call interp_met_field(input_name, fg_data%field, HH, M, &
1356                                      field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1357                                      we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1358                                      halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
1359                                      field%modified_mask, process_bdy_width, landmask)
1360                      end if
1361                   else
1362                      call mprintf(.true.,INFORM,' - already processed this field from constant file.')
1363                   end if
1364
1365                end if
1366
1367                call bitarray_merge(field%valid_mask, field%modified_mask)
1368
1369                deallocate(halo_slab)
1370
1371                ! Store the interpolated field
1372                call storage_put_field(field)
1373
1374                call pop_source_projection()
1375
1376             end if
1377          end do
1378
1379          call read_met_close()
1380
1381          call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, &
1382                                      fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, &
1383                                      fg_data%deltalon, fg_data%starti, fg_data%startj, &
1384                                      fg_data%startlat, fg_data%startlon,  &
1385                                      fg_data%pole_lat, fg_data%pole_lon,        &
1386                                      fg_data%centerlat, fg_data%centerlon,      &
1387                                      real(fg_data%nx+1)/2., real(fg_data%ny+1)/2.,  &
1388                                      earth_radius=fg_data%earth_radius*1000.)
1389
1390          !
1391          ! If necessary, rotate winds to earth-relative for this fg source
1392          !
1393
1394          call storage_get_levels(u_field, u_levels)
1395          call storage_get_levels(v_field, v_levels)
1396
1397          if (associated(u_levels) .and. associated(v_levels)) then
1398             u_idx = 1
1399             do u_idx = 1, size(u_levels)
1400                u_field%header%vertical_level = u_levels(u_idx)
1401                call storage_get_field(u_field, istatus)
1402                v_field%header%vertical_level = v_levels(u_idx)
1403                call storage_get_field(v_field, istatus)
1404
1405                if (associated(u_field%modified_mask) .and. &
1406                    associated(v_field%modified_mask)) then
1407
1408                   if (u_field%header%is_wind_grid_rel) then
1409                      if (gridtype == 'C') then
1410                         call map_to_met(u_field%r_arr, u_field%modified_mask, &
1411                                         v_field%r_arr, v_field%modified_mask, &
1412                                         we_mem_stag_s, sn_mem_s, &
1413                                         we_mem_stag_e, sn_mem_e, &
1414                                         we_mem_s, sn_mem_stag_s, &
1415                                         we_mem_e, sn_mem_stag_e, &
1416                                         xlon_u, xlon_v, xlat_u, xlat_v)
1417                      else if (gridtype == 'E') then
1418                         call map_to_met_nmm(u_field%r_arr, u_field%modified_mask, &
1419                                             v_field%r_arr, v_field%modified_mask, &
1420                                             we_mem_s, sn_mem_s, &
1421                                             we_mem_e, sn_mem_e, &
1422                                             xlat_v, xlon_v)
1423                      end if
1424                   end if
1425
1426                   call bitarray_destroy(u_field%modified_mask)
1427                   call bitarray_destroy(v_field%modified_mask)
1428                   nullify(u_field%modified_mask)
1429                   nullify(v_field%modified_mask)
1430                   call storage_put_field(u_field)
1431                   call storage_put_field(v_field)
1432                end if
1433
1434             end do
1435
1436             deallocate(u_levels)
1437             deallocate(v_levels)
1438
1439          end if
1440
1441          call pop_source_projection()
1442
1443       else
1444          if (do_const_processing) then
1445             call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=input_name)
1446          else
1447             call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=trim(input_name)//':'//trim(temp_date))
1448          end if
1449       end if
1450
1451    end subroutine process_intermediate_fields
1452
1453
1454    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1455    ! Name: process_mpas_fields
1456    !
1457    ! Purpose:
1458    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1459    subroutine process_mpas_fields(input_name, do_const_processing, temp_date, fg_data, got_this_field, &
1460                                           landmask, process_bdy_width, &
1461                                           u_field, v_field, &
1462                                           dom_dx, dom_dy, &
1463                                           xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v, &
1464                                           output_flags, west_east_dim, south_north_dim, &
1465                                           we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1466                                           we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
1467                                           sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
1468                                           we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
1469                                           sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e )
1470
1471       use bitarray_module
1472       use gridinfo_module
1473       use interp_module
1474       use interp_option_module
1475       use read_met_module
1476       use llxy_module
1477       use misc_definitions_module
1478       use module_debug
1479       use output_module
1480       use parallel_module
1481       use rotate_winds_module
1482       use storage_module
1483       use scan_input
1484       use mpas_mesh
1485
1486       implicit none
1487
1488 ! BUG: Move this constant to misc_definitions_module?
1489 integer, parameter :: BDR_WIDTH = 3
1490
1491       character (len=*), intent(inout) :: input_name
1492       logical, intent(in) :: do_const_processing
1493       character (len=*), intent(in) :: temp_date
1494       type (met_data), intent(inout) :: fg_data
1495       logical, dimension(:), intent(inout) :: got_this_field
1496       real, pointer, dimension(:,:) :: landmask
1497       integer, intent(in) :: process_bdy_width
1498       type (fg_input), intent(inout) :: u_field, v_field
1499       character (len=128), dimension(:), intent(inout) :: output_flags
1500       real, intent(in) :: dom_dx, dom_dy
1501       real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
1502       integer, intent(in) :: west_east_dim, south_north_dim, &
1503                  we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1504                  we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
1505                  sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
1506                  we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
1507                  sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e
1508
1509       real, parameter :: deg2rad = asin(1.0) / 90.0
1510
1511       integer :: i, j, k
1512       integer :: idx
1513       integer :: istat
1514       integer :: strlen
1515       character (len=MAX_FILENAME_LEN) :: mpas_filename
1516       integer :: nRecords
1517       type (input_handle_type) :: mpas_handle
1518       type (input_field_type) :: mpas_field
1519       type (target_field_type) :: wrf_field
1520       type (fg_input) :: field_to_store
1521
1522       strlen = len_trim(input_name)
1523       if (do_const_processing) then
1524          write(mpas_filename,'(a)') input_name(6:strlen)
1525       else
1526          write(mpas_filename,'(a)') input_name(6:strlen)//'.'//trim(temp_date)//'.nc'
1527       end if
1528       call mprintf(.true.,LOGFILE,'Processing MPAS fields from file %s',s1=mpas_filename)
1529
1530       !
1531       ! If we do not already have mesh information, get that now...
1532       !
1533       if (.not. mpas_source_mesh % valid) then
1534          if (mpas_mesh_setup(mpas_filename, mpas_source_mesh) /= 0) then
1535             call mprintf(.true.,ERROR, 'Error setting up MPAS mesh %s with scan_input_open', s1=mpas_filename)
1536          end if
1537       end if
1538
1539       !
1540       ! If we have not already defined the WRF grid, do that now...
1541       !
1542       if (.not. wrf_target_mesh_m % valid) then
1543          allocate(xlat_rad(size(xlat,1), size(xlat,2)))
1544          allocate(xlon_rad(size(xlat,1), size(xlat,2)))
1545          xlat_rad(:,:) = xlat(:,:) * deg2rad
1546          xlon_rad(:,:) = xlon(:,:) * deg2rad
1547          call mprintf(.true.,LOGFILE,'Need to set up WRF target mass-grid')
1548          if (target_mesh_setup(wrf_target_mesh_m, lat2d=xlat_rad, lon2d=xlon_rad) /= 0) then
1549             call mprintf(.true.,ERROR, 'Error setting up WRF target grid')
1550          end if
1551
1552          call mprintf(.true.,LOGFILE,'Also computing remapping weights...')
1553          if (remap_info_setup(mpas_source_mesh, wrf_target_mesh_m, remap_info_m) /= 0) then
1554             call mprintf(.true.,ERROR, 'Error computing remapping weights from MPAS to WRF grid')
1555          end if
1556       else
1557          call mprintf(.true.,LOGFILE,'Already set up WRF target mass-grid')
1558       end if
1559
1560       if (.not. wrf_target_mesh_u % valid) then
1561          allocate(xlat_u_rad(size(xlat_u,1), size(xlat_u,2)))
1562          allocate(xlon_u_rad(size(xlat_u,1), size(xlat_u,2)))
1563          xlat_u_rad(:,:) = xlat_u(:,:) * deg2rad
1564          xlon_u_rad(:,:) = xlon_u(:,:) * deg2rad
1565          call mprintf(.true.,LOGFILE,'Need to set up WRF target U-grid')
1566          if (target_mesh_setup(wrf_target_mesh_u, lat2d=xlat_u_rad, lon2d=xlon_u_rad) /= 0) then
1567             call mprintf(.true.,ERROR, 'Error setting up WRF target grid')
1568          end if
1569
1570          call mprintf(.true.,LOGFILE,'Also computing remapping weights...')
1571          if (remap_info_setup(mpas_source_mesh, wrf_target_mesh_u, remap_info_u) /= 0) then
1572             call mprintf(.true.,ERROR, 'Error computing remapping weights from MPAS to WRF grid')
1573          end if
1574       else
1575          call mprintf(.true.,LOGFILE,'Already set up WRF target U-grid')
1576       end if
1577
1578       if (.not. wrf_target_mesh_v % valid) then
1579          allocate(xlat_v_rad(size(xlat_v,1), size(xlat_v,2)))
1580          allocate(xlon_v_rad(size(xlat_v,1), size(xlat_v,2)))
1581          xlat_v_rad(:,:) = xlat_v(:,:) * deg2rad
1582          xlon_v_rad(:,:) = xlon_v(:,:) * deg2rad
1583          call mprintf(.true.,LOGFILE,'Need to set up WRF target V-grid')
1584          if (target_mesh_setup(wrf_target_mesh_v, lat2d=xlat_v_rad, lon2d=xlon_v_rad) /= 0) then
1585             call mprintf(.true.,ERROR, 'Error setting up WRF target grid')
1586          end if
1587
1588          call mprintf(.true.,LOGFILE,'Also computing remapping weights...')
1589          if (remap_info_setup(mpas_source_mesh, wrf_target_mesh_v, remap_info_v) /= 0) then
1590             call mprintf(.true.,ERROR, 'Error computing remapping weights from MPAS to WRF grid')
1591          end if
1592       else
1593          call mprintf(.true.,LOGFILE,'Already set up WRF target V-grid')
1594       end if
1595
1596
1597       if (scan_input_open(mpas_filename, mpas_handle, nRecords) /= 0) then
1598          call mprintf(.true.,ERROR, 'Error opening %s with scan_input_open', s1=mpas_filename)
1599       end if
1600
1601
1602       ! Initialize fg_input structure to store the field
1603       field_to_store%header%version = 1
1604       field_to_store%header%date = '?'
1605       if (do_const_processing) then
1606          field_to_store%header%time_dependent = .false.
1607          field_to_store%header%constant_field = .true.
1608       else
1609          field_to_store%header%time_dependent = .true.
1610          field_to_store%header%constant_field = .false.
1611       end if
1612       field_to_store%header%forecast_hour = 0.0
1613       field_to_store%header%fg_source = 'FG'
1614       field_to_store%header%field = ' '
1615       field_to_store%header%field(1:9) = '?'
1616       field_to_store%header%units = ' '
1617       field_to_store%header%units(1:25) = '?'
1618       field_to_store%header%description = ' '
1619       field_to_store%header%description(1:46) = '?'
1620       field_to_store%header%vertical_coord = 'z_dim_name'
1621       field_to_store%header%vertical_level = 0
1622       field_to_store%header%sr_x = 1
1623       field_to_store%header%sr_y = 1
1624       field_to_store%header%array_order = 'XY '
1625       field_to_store%header%is_wind_grid_rel = .false.
1626       field_to_store%header%array_has_missing_values = .false.
1627       nullify(field_to_store%r_arr)
1628       nullify(field_to_store%valid_mask)
1629       nullify(field_to_store%modified_mask)
1630
1631       ! If we should not output this field, just list it as a mask field
1632 !???      if (output_this_field(idx)) then
1633          field_to_store%header%mask_field = .false.
1634 !???      else
1635 !???         field%header%mask_field = .true.
1636 !???      end if
1637
1638
1639       do while (scan_input_next_field(mpas_handle, mpas_field) == 0)
1640
1641          if (can_remap_field(mpas_field)) then
1642
1643             ! Here, rename a few MPAS fields that would be difficult to treat
1644             ! with METGRID.TBL options; principally, these are surface fields
1645             ! that have different names from their upper-air counterparts.
1646             if (trim(mpas_field % name) == 'u10') then
1647                mpas_field % name = 'uReconstructZonal'
1648             else if (trim(mpas_field % name) == 'v10') then
1649                mpas_field % name = 'uReconstructMeridional'
1650             else if (trim(mpas_field % name) == 'q2') then
1651                mpas_field % name = 'qv'
1652             else if (trim(mpas_field % name) == 't2m') then
1653                mpas_field % name = 'theta'
1654             end if
1655
1656             ! Mark this MPAS field as "gotten" for any later checks
1657             ! on mandatory fields
1658             idx = mpas_name_to_idx(trim(mpas_field % name))
1659             if (idx > 0) then
1660                got_this_field(idx) = .true.
1661                if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
1662                   output_flags(idx) = flag_in_output(idx)
1663                end if
1664             end if
1665
1666             istat = scan_input_read_field(mpas_field, frame=1)
1667
1668             field_to_store%map%stagger = mpas_output_stagger(mpas_field % name)
1669             if (field_to_store%map%stagger == M) then
1670                field_to_store%header%dim1(1) = we_mem_s
1671                field_to_store%header%dim1(2) = we_mem_e
1672                field_to_store%header%dim2(1) = sn_mem_s
1673                field_to_store%header%dim2(2) = sn_mem_e
1674                if (idx > 0) then
1675                   if (masked(idx) == MASKED_WATER) then
1676                      istat = remap_field(remap_info_m, mpas_field, wrf_field, masked=.true.)
1677                   else
1678                      istat = remap_field(remap_info_m, mpas_field, wrf_field, masked=.false.)
1679                   end if
1680                else
1681                   istat = remap_field(remap_info_m, mpas_field, wrf_field, masked=.false.)
1682                end if
1683             else if (field_to_store%map%stagger == U) then
1684                field_to_store%header%dim1(1) = we_mem_stag_s
1685                field_to_store%header%dim1(2) = we_mem_stag_e
1686                field_to_store%header%dim2(1) = sn_mem_s
1687                field_to_store%header%dim2(2) = sn_mem_e
1688                istat = remap_field(remap_info_u, mpas_field, wrf_field)
1689             else if (field_to_store%map%stagger == V) then
1690                field_to_store%header%dim1(1) = we_mem_s
1691                field_to_store%header%dim1(2) = we_mem_e
1692                field_to_store%header%dim2(1) = sn_mem_stag_s
1693                field_to_store%header%dim2(2) = sn_mem_stag_e
1694                istat = remap_field(remap_info_v, mpas_field, wrf_field)
1695             else
1696                call mprintf(.true.,ERROR, 'Cannot handle requested output stagger %i for MPAS field %s ...', &
1697                             i1=field_to_store%map%stagger, s1=trim(mpas_field % name))
1698             end if
1699
1700             if (wrf_field % ndims == 3 .and. trim(wrf_field % dimnames(1)) == 'nVertLevels') then   ! 3-d MPAS atmosphere field
1701                field_to_store % header % field = mpas_to_intermediate_name(mpas_field % name)
1702
1703                ! If no match in the METGRID.TBL was found for this MPAS field, just use the MPAS name
1704                if (len_trim(field_to_store % header % field) == 0) then
1705                   field_to_store % header % field = trim(mpas_field % name)
1706                end if
1707
1708                field_to_store % header % vertical_coord = 'num_mpas_levels'
1709                do k=1,wrf_field % dimlens(1)
1710                   allocate(field_to_store % r_arr(wrf_field % dimlens(2), wrf_field % dimlens(3)))
1711                   allocate(field_to_store % valid_mask)
1712                   call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(2), wrf_field % dimlens(3))
1713                   do j=1,wrf_field % dimlens(3)
1714                   do i=1,wrf_field % dimlens(2)
1715                      call bitarray_set(field_to_store % valid_mask, i, j)
1716                   end do
1717                   end do
1718                   field_to_store % header % vertical_level = k
1719
1720                   if (wrf_field % xtype == FIELD_TYPE_REAL) then
1721                      field_to_store % r_arr(:,:) = wrf_field % array3r(k,:,:)
1722                   else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then
1723                      field_to_store % r_arr(:,:) = wrf_field % array3d(k,:,:)
1724                   end if
1725
1726                   ! The u_field and v_field are used later by calling code to
1727                   ! determine which fields represent the x- and y-components of
1728                   ! horizonal wind velocity for the purposes of wind rotation
1729                   if (trim(mpas_field % name) == 'uReconstructZonal') then
1730                      call dup(field_to_store, u_field)
1731                   end if
1732                   if (trim(mpas_field % name) == 'uReconstructMeridional') then
1733                      call dup(field_to_store, v_field)
1734                   end if
1735
1736                   call storage_put_field(field_to_store)
1737                end do
1738
1739             else if (wrf_field % ndims == 3 .and. trim(wrf_field % dimnames(1)) == 'nVertLevelsP1') then   ! 3-d MPAS atmosphere field
1740                field_to_store % header % field = mpas_to_intermediate_name(mpas_field % name)
1741
1742                ! If no match in the METGRID.TBL was found for this MPAS field, just use the MPAS name
1743                if (len_trim(field_to_store % header % field) == 0) then
1744                   field_to_store % header % field = trim(mpas_field % name)
1745                end if
1746
1747                ! Handle surface level
1748                field_to_store % header % vertical_coord = 'none'
1749                allocate(field_to_store % r_arr(wrf_field % dimlens(2), wrf_field % dimlens(3)))
1750                allocate(field_to_store % valid_mask)
1751                call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(2), wrf_field % dimlens(3))
1752                do j=1,wrf_field % dimlens(3)
1753                do i=1,wrf_field % dimlens(2)
1754                   call bitarray_set(field_to_store % valid_mask, i, j)
1755                end do
1756                end do
1757                field_to_store % header % vertical_level = 200100.0
1758
1759                if (wrf_field % xtype == FIELD_TYPE_REAL) then
1760                   field_to_store % r_arr(:,:) = wrf_field % array3r(1,:,:)
1761                else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then
1762                   field_to_store % r_arr(:,:) = wrf_field % array3d(1,:,:)
1763                end if
1764
1765                call storage_put_field(field_to_store)
1766
1767                ! Handle all other layers
1768                field_to_store % header % vertical_coord = 'num_mpas_levels'
1769                do k=1,wrf_field % dimlens(1) - 1
1770                   allocate(field_to_store % r_arr(wrf_field % dimlens(2), wrf_field % dimlens(3)))
1771                   allocate(field_to_store % valid_mask)
1772                   call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(2), wrf_field % dimlens(3))
1773                   do j=1,wrf_field % dimlens(3)
1774                   do i=1,wrf_field % dimlens(2)
1775                      call bitarray_set(field_to_store % valid_mask, i, j)
1776                   end do
1777                   end do
1778                   field_to_store % header % vertical_level = k
1779
1780                   ! Average to layer midpoint
1781                   if (wrf_field % xtype == FIELD_TYPE_REAL) then
1782                      field_to_store % r_arr(:,:) = 0.5 * (wrf_field % array3r(k,:,:) + wrf_field % array3r(k+1,:,:))
1783                   else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then
1784                      field_to_store % r_arr(:,:) = 0.5 * (wrf_field % array3d(k,:,:) + wrf_field % array3d(k+1,:,:))
1785                   end if
1786
1787                   call storage_put_field(field_to_store)
1788                end do
1789
1790                ! Special case: zgrid field also provides SOILHGT
1791                if (trim(mpas_field % name) == 'zgrid') then
1792                   field_to_store % header % field = 'SOILHGT'
1793
1794                   allocate(field_to_store % r_arr(wrf_field % dimlens(2), wrf_field % dimlens(3)))
1795                   allocate(field_to_store % valid_mask)
1796                   call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(2), wrf_field % dimlens(3))
1797                   do j=1,wrf_field % dimlens(3)
1798                   do i=1,wrf_field % dimlens(2)
1799                      call bitarray_set(field_to_store % valid_mask, i, j)
1800                   end do
1801                   end do
1802                   field_to_store % header % vertical_level = 200100.0
1803
1804                   if (wrf_field % xtype == FIELD_TYPE_REAL) then
1805                      field_to_store % r_arr(:,:) = wrf_field % array3r(1,:,:)
1806                   else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then
1807                      field_to_store % r_arr(:,:) = wrf_field % array3d(1,:,:)
1808                   end if
1809
1810                   call storage_put_field(field_to_store)
1811
1812                   do idx=1,num_entries
1813                      if (trim(fieldname(idx)) == 'SOILHGT') then
1814                         got_this_field(idx) = .true.
1815                         if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
1816                            output_flags(idx) = flag_in_output(idx)
1817                         end if
1818                         exit
1819                      end if
1820                   end do
1821                end if
1822
1823             else if (wrf_field % ndims == 3 .and. trim(wrf_field % dimnames(1)) == 'nSoilLevels') then   ! 3-d MPAS soil field
1824
1825                field_to_store % header % vertical_coord = 'none'
1826                if (trim(mpas_field % name) == 'tslb') then
1827                   do k=1,wrf_field % dimlens(1)
1828                      if (k == 1) then
1829                         field_to_store % header % field = 'ST000010'
1830                      else if (k == 2) then
1831                         field_to_store % header % field = 'ST010040'
1832                      else if (k == 3) then
1833                         field_to_store % header % field = 'ST040100'
1834                      else if (k == 4) then
1835                         field_to_store % header % field = 'ST100200'
1836                      else
1837                         call mprintf(.true.,ERROR, 'Too many soil layers in MPAS soil field %s ...', s1=trim(mpas_field % name))
1838                      end if
1839                      allocate(field_to_store % r_arr(wrf_field % dimlens(2), wrf_field % dimlens(3)))
1840                      allocate(field_to_store % valid_mask)
1841                      call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(2), wrf_field % dimlens(3))
1842                      do j=1,wrf_field % dimlens(3)
1843                      do i=1,wrf_field % dimlens(2)
1844                         call bitarray_set(field_to_store % valid_mask, i, j)
1845                      end do
1846                      end do
1847                      field_to_store % header % vertical_level = 200100.0
1848
1849                      if (wrf_field % xtype == FIELD_TYPE_REAL) then
1850                         field_to_store % r_arr(:,:) = wrf_field % array3r(k,:,:)
1851                      else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then
1852                         field_to_store % r_arr(:,:) = wrf_field % array3d(k,:,:)
1853                      end if
1854
1855                      if (idx > 0) then
1856                         if (masked(idx) == MASKED_WATER) then
1857                            where (landmask(:,:) == 0) field_to_store % r_arr(:,:) = fill_missing(idx)
1858                         end if
1859                      end if
1860
1861                      call storage_put_field(field_to_store)
1862                   end do
1863                else if (trim(mpas_field % name) == 'smois') then
1864                   do k=1,wrf_field % dimlens(1)
1865                      if (k == 1) then
1866                         field_to_store % header % field = 'SM000010'
1867                      else if (k == 2) then
1868                         field_to_store % header % field = 'SM010040'
1869                      else if (k == 3) then
1870                         field_to_store % header % field = 'SM040100'
1871                      else if (k == 4) then
1872                         field_to_store % header % field = 'SM100200'
1873                      else
1874                         call mprintf(.true.,ERROR, 'Too many soil layers in MPAS soil field %s ...', s1=trim(mpas_field % name))
1875                      end if
1876                      allocate(field_to_store % r_arr(wrf_field % dimlens(2), wrf_field % dimlens(3)))
1877                      allocate(field_to_store % valid_mask)
1878                      call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(2), wrf_field % dimlens(3))
1879                      do j=1,wrf_field % dimlens(3)
1880                      do i=1,wrf_field % dimlens(2)
1881                         call bitarray_set(field_to_store % valid_mask, i, j)
1882                      end do
1883                      end do
1884                      field_to_store % header % vertical_level = 200100.0
1885
1886                      if (wrf_field % xtype == FIELD_TYPE_REAL) then
1887                         field_to_store % r_arr(:,:) = wrf_field % array3r(k,:,:)
1888                      else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then
1889                         field_to_store % r_arr(:,:) = wrf_field % array3d(k,:,:)
1890                      end if
1891
1892                      if (idx > 0) then
1893                         if (masked(idx) == MASKED_WATER) then
1894                            where (landmask(:,:) == 0) field_to_store % r_arr(:,:) = fill_missing(idx)
1895                         end if
1896                      end if
1897
1898                      call storage_put_field(field_to_store)
1899                   end do
1900                else
1901                   call mprintf(.true.,WARN, 'Skipping unknown MPAS soil field %s ...', s1=trim(mpas_field % name))
1902                end if
1903
1904             else if (wrf_field % ndims == 2) then   ! 2-d MPAS field
1905                field_to_store % header % field = mpas_to_intermediate_name(mpas_field % name)
1906
1907                ! If no match in the METGRID.TBL was found for this MPAS field, just use the MPAS name
1908                if (len_trim(field_to_store % header % field) == 0) then
1909                   field_to_store % header % field = trim(mpas_field % name)
1910                end if
1911
1912                field_to_store % header % vertical_coord = 'none'
1913                allocate(field_to_store % r_arr(wrf_field % dimlens(1), wrf_field % dimlens(2)))
1914                allocate(field_to_store % valid_mask)
1915                call bitarray_create(field_to_store % valid_mask, wrf_field % dimlens(1), wrf_field % dimlens(2))
1916                do j=1,wrf_field % dimlens(2)
1917                do i=1,wrf_field % dimlens(1)
1918                   call bitarray_set(field_to_store % valid_mask, i, j)
1919                end do
1920                end do
1921                field_to_store % header % vertical_level = 200100.0
1922
1923                if (wrf_field % xtype == FIELD_TYPE_REAL) then
1924                   field_to_store % r_arr(:,:) = wrf_field % array2r(:,:)
1925                else if (wrf_field % xtype == FIELD_TYPE_DOUBLE) then
1926                   field_to_store % r_arr(:,:) = wrf_field % array2d(:,:)
1927                end if
1928
1929                if (idx > 0) then
1930                   if (masked(idx) == MASKED_WATER) then
1931                      where (landmask(:,:) == 0) field_to_store % r_arr(:,:) = fill_missing(idx)
1932                   end if
1933                end if
1934
1935                call storage_put_field(field_to_store)
1936             end if
1937
1938             istat = free_target_field(wrf_field)
1939          end if
1940
1941          istat = scan_input_free_field(mpas_field)
1942       end do
1943
1944       if (scan_input_close(mpas_handle) /= 0) then
1945          call mprintf(.true.,ERROR, 'Error closing %s with scan_input_close', s1=mpas_filename)
1946       end if
1947
1948    end subroutine process_mpas_fields
1949
1950
1951    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1952    ! Name: derive_mpas_fields
1953    !
1954    ! Purpose:
1955    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1956    subroutine derive_mpas_fields()
1957
1958       use bitarray_module
1959       use gridinfo_module
1960       use interp_module
1961       use interp_option_module
1962       use read_met_module
1963       use llxy_module
1964       use misc_definitions_module
1965       use module_debug
1966       use output_module
1967       use parallel_module
1968       use rotate_winds_module
1969       use storage_module
1970       use scan_input
1971       use mpas_mesh
1972       use constants_module
1973
1974       implicit none
1975
1976       integer :: k
1977       integer :: istatus
1978       integer, pointer, dimension(:) :: theta_levels => null()
1979       integer, pointer, dimension(:) :: pressure_levels => null()
1980       real, pointer, dimension(:,:) :: exner => null()
1981       type (fg_input) :: theta_field
1982       type (fg_input) :: pressure_field
1983
1984       !
1985       ! Derive TT from theta and pressure
1986       !
1987       theta_field%header%time_dependent = .true.
1988       theta_field%header%constant_field = .false.
1989       theta_field%header%field = 'TT'
1990       theta_field%header%vertical_level = 0
1991       nullify(theta_field%r_arr)
1992       nullify(theta_field%valid_mask)
1993       nullify(theta_field%modified_mask)
1994
1995       pressure_field%header%time_dependent = .true.
1996       pressure_field%header%constant_field = .false.
1997       pressure_field%header%field = 'PRESSURE'
1998       pressure_field%header%vertical_level = 0
1999       nullify(pressure_field%r_arr)
2000       nullify(pressure_field%valid_mask)
2001       nullify(pressure_field%modified_mask)
2002
2003       call storage_get_levels(theta_field, theta_levels)
2004       call storage_get_levels(pressure_field, pressure_levels)
2005
2006       if (associated(theta_levels) .and. associated(pressure_levels)) then
2007          call mprintf(.true.,LOGFILE, 'Computing MPAS TT field from theta and pressure...')
2008
2009          if (size(theta_levels) == size(pressure_levels)) then
2010             do k = 1, size(theta_levels)
2011                call mprintf(.true.,LOGFILE, 'Computing TT at level %i for MPAS dataset', i1=theta_levels(k))
2012                theta_field % header % vertical_level = theta_levels(k)
2013                call storage_get_field(theta_field, istatus)
2014                if (istatus /= 0) then
2015                   call mprintf(.true.,ERROR, 'Could not get MPAS theta field at level %i', i1=theta_levels(k))
2016                   return
2017                end if
2018
2019                pressure_field % header % vertical_level = pressure_levels(k)
2020                call storage_get_field(pressure_field, istatus)
2021                if (istatus /= 0) then
2022                   call mprintf(.true.,ERROR, 'Could not get MPAS pressure field at level %i', i1=theta_levels(k))
2023                   return
2024                end if
2025
2026                ! Compute temperature
2027                if (.not. associated(exner)) then
2028                   allocate(exner(size(theta_field % r_arr, 1), size(theta_field % r_arr, 2)))
2029                end if
2030                exner(:,:) = (pressure_field % r_arr(:,:) / P0)**(RD/CP)
2031                theta_field % r_arr(:,:) = theta_field % r_arr(:,:) * exner(:,:)
2032
2033                call storage_put_field(theta_field)
2034             end do
2035          else
2036             call mprintf(.true.,ERROR, 'The MPAS theta and pressure fields do not have the same number of levels!')
2037          end if
2038
2039          deallocate(theta_levels)
2040          deallocate(pressure_levels)
2041          if (associated(exner)) then
2042             deallocate(exner)
2043          end if
2044       end if
2045
2046    end subroutine derive_mpas_fields
2047
2048
2049    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2050    ! Name: get_interp_masks
2051    !
2052    ! Purpose:
2053    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2054    subroutine get_interp_masks(fg_prefix, is_constants, fg_date)
2055
2056       use interp_option_module
2057       use read_met_module
2058       use storage_module
2059
2060       implicit none
2061
2062       ! Arguments
2063       logical, intent(in) :: is_constants
2064       character (len=*), intent(in) :: fg_prefix, fg_date
2065
2066 ! BUG: Move this constant to misc_definitions_module?
2067 integer, parameter :: BDR_WIDTH = 3
2068
2069       ! Local variables
2070       integer :: i, istatus, idx, idxt
2071       type (fg_input) :: mask_field
2072       type (met_data) :: fg_data
2073
2074       istatus = 0
2075
2076       call read_met_init(fg_prefix, is_constants, fg_date, istatus)
2077
2078       do while (istatus == 0)
2079
2080          call read_next_met_field(fg_data, istatus)
2081
2082          if (istatus == 0) then
2083
2084             ! Find out which METGRID.TBL entry goes with this field
2085             idxt = num_entries + 1
2086             do idx=1,num_entries
2087                if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
2088                    (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
2089
2090                   if (index(fg_prefix,trim(from_input(idx))) /= 0 .or. &
2091                      (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
2092                      idxt = idx
2093                   end if
2094
2095                end if
2096             end do
2097             idx = idxt
2098             if (idx > num_entries) idx = num_entries ! The last entry is a default
2099
2100             ! Do we need to rename this field?
2101             if (output_name(idx) /= ' ') then
2102                fg_data%field = output_name(idx)(1:9)
2103
2104                idxt = num_entries + 1
2105                do idx=1,num_entries
2106                   if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
2107                       (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
2108
2109                      if (index(fg_prefix,trim(from_input(idx))) /= 0 .or. &
2110                         (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
2111                         idxt = idx
2112                      end if
2113
2114                   end if
2115                end do
2116                idx = idxt
2117                if (idx > num_entries) idx = num_entries ! The last entry is a default
2118             end if
2119
2120             do i=1,num_entries
2121                if (interp_mask(i) /= ' ' .and. (trim(interp_mask(i)) == trim(fg_data%field))) then
2122
2123                   mask_field%header%version = 1
2124                   mask_field%header%date = ' '
2125                   mask_field%header%date = fg_date
2126                   if (is_constants) then
2127                      mask_field%header%time_dependent = .false.
2128                      mask_field%header%constant_field = .true.
2129                   else
2130                      mask_field%header%time_dependent = .true.
2131                      mask_field%header%constant_field = .false.
2132                   end if
2133                   mask_field%header%mask_field = .true.
2134                   mask_field%header%forecast_hour = 0.
2135                   mask_field%header%fg_source = 'degribbed met data'
2136                   mask_field%header%field = trim(fg_data%field)//'.mask'
2137                   mask_field%header%units = '-'
2138                   mask_field%header%description = '-'
2139                   mask_field%header%vertical_coord = 'none'
2140                   mask_field%header%vertical_level = 1
2141                   mask_field%header%sr_x = 1
2142                   mask_field%header%sr_y = 1
2143                   mask_field%header%array_order = 'XY'
2144                   mask_field%header%dim1(1) = 1
2145                   mask_field%header%dim1(2) = fg_data%nx
2146                   mask_field%header%dim2(1) = 1
2147                   mask_field%header%dim2(2) = fg_data%ny
2148                   mask_field%header%is_wind_grid_rel = .true.
2149                   mask_field%header%array_has_missing_values = .false.
2150                   mask_field%map%stagger = M
2151
2152                   ! Do a simple check to see whether this is a global lat/lon dataset
2153                   ! Note that we do not currently support regional Gaussian grids
2154                   if ((fg_data%iproj == PROJ_LATLON .and. abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) &
2155                        .or. (fg_data%iproj == PROJ_GAUSS)) then
2156                      allocate(mask_field%r_arr(1-BDR_WIDTH:fg_data%nx+BDR_WIDTH,1:fg_data%ny))
2157
2158                      mask_field%r_arr(1:fg_data%nx,                      1:fg_data%ny) = &
2159                          fg_data%slab(1:fg_data%nx,              1:fg_data%ny)
2160
2161                      mask_field%r_arr(1-BDR_WIDTH:0,                     1:fg_data%ny) = &
2162                          fg_data%slab(fg_data%nx-BDR_WIDTH+1:fg_data%nx, 1:fg_data%ny)
2163
2164                      mask_field%r_arr(fg_data%nx+1:fg_data%nx+BDR_WIDTH, 1:fg_data%ny) = &
2165                          fg_data%slab(1:BDR_WIDTH,       1:fg_data%ny)
2166
2167                   else
2168                      allocate(mask_field%r_arr(1:fg_data%nx,1:fg_data%ny))
2169                      mask_field%r_arr = fg_data%slab
2170                   end if
2171
2172                   nullify(mask_field%valid_mask)
2173                   nullify(mask_field%modified_mask)
2174
2175                   call storage_put_field(mask_field)
2176
2177                   exit
2178
2179                end if
2180             end do
2181
2182             if (associated(fg_data%slab)) deallocate(fg_data%slab)
2183
2184          end if
2185
2186       end do
2187
2188       call read_met_close()
2189
2190    end subroutine get_interp_masks
2191
2192
2193    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2194    ! Name: interp_met_field
2195    !
2196    ! Purpose:
2197    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2198    subroutine interp_met_field(input_name, short_fieldnm, ifieldstagger, istagger, &
2199                                field, xlat, xlon, sm1, em1, sm2, em2, &
2200                                sd1, ed1, sd2, ed2, &
2201                                slab, minx, maxx, miny, maxy, bdr, do_gcell_interp, &
2202                                new_pts, process_bdy_width, landmask)
2203
2204       use bitarray_module
2205       use interp_module
2206       use interp_option_module
2207       use llxy_module
2208       use misc_definitions_module
2209       use storage_module
2210
2211       implicit none
2212
2213       ! Arguments
2214       integer, intent(in) :: ifieldstagger, istagger, &
2215                              sm1, em1, sm2, em2, &
2216                              sd1, ed1, sd2, ed2, &
2217                              minx, maxx, miny, maxy, bdr, &
2218                              process_bdy_width
2219       real, dimension(minx:maxx,miny:maxy), intent(in) :: slab
2220       real, dimension(sm1:em1,sm2:em2), intent(in) :: xlat, xlon
2221       real, dimension(sm1:em1,sm2:em2), intent(in), optional :: landmask
2222       logical, intent(in) :: do_gcell_interp
2223       character (len=9), intent(in) :: short_fieldnm
2224       character (len=MAX_FILENAME_LEN), intent(in) :: input_name
2225       type (fg_input), intent(inout) :: field
2226       type (bitarray), intent(inout) :: new_pts
2227
2228       ! Local variables
2229       integer :: i, j, idx, idxt, orig_selected_proj, interp_mask_status, &
2230                  interp_land_mask_status, interp_water_mask_status, process_width
2231       integer, pointer, dimension(:) :: interp_array, interp_opts
2232       real :: rx, ry, temp
2233       real, pointer, dimension(:,:) :: data_count
2234       type (fg_input) :: mask_field, mask_water_field, mask_land_field
2235       !BPR BEGIN
2236       real, dimension(sm1:em1,sm2:em2) :: r_arr_cur_source
2237       !BPR END
2238
2239       ! CWH Initialize local pointer variables
2240       nullify(interp_array)
2241       nullify(interp_opts)
2242       nullify(data_count)
2243
2244       ! Find index into fieldname, interp_method, masked, and fill_missing
2245       !   of the current field
2246       idxt = num_entries + 1
2247       do idx=1,num_entries
2248          if ((index(fieldname(idx), trim(short_fieldnm)) /= 0) .and. &
2249              (len_trim(fieldname(idx)) == len_trim(short_fieldnm))) then
2250             if (index(input_name,trim(from_input(idx))) /= 0 .or. &
2251                (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
2252                idxt = idx
2253             end if
2254          end if
2255       end do
2256       idx = idxt
2257       if (idx > num_entries) then
2258          call mprintf(.true.,WARN,'Entry in METGRID.TBL not found for field %s. '// &
2259                       'Default options will be used for this field!', s1=short_fieldnm)
2260          idx = num_entries ! The last entry is a default
2261       end if
2262
2263       if (process_bdy_width == 0) then
2264          process_width = max(ed1-sd1+1, ed2-sd2+1)
2265       else
2266          process_width = process_bdy_width
2267         ! Add two extra rows/cols to accommodate staggered fields: one extra row/col for
2268         !    averaging to mass points in real, and one beyond that for averaging during
2269         !    wind rotation
2270         if (ifieldstagger /= M) process_width = process_width + 2
2271       end if
2272
2273       field%header%dim1(1) = sm1
2274       field%header%dim1(2) = em1
2275       field%header%dim2(1) = sm2
2276       field%header%dim2(2) = em2
2277       field%map%stagger = ifieldstagger
2278       if (.not. associated(field%r_arr)) then
2279          allocate(field%r_arr(sm1:em1,sm2:em2))
2280       end if
2281
2282       interp_mask_status = 1
2283       interp_land_mask_status = 1
2284       interp_water_mask_status = 1
2285
2286       if (interp_mask(idx) /= ' ') then
2287          mask_field%header%version = 1
2288          mask_field%header%forecast_hour = 0.
2289          mask_field%header%field = trim(interp_mask(idx))//'.mask'
2290          mask_field%header%vertical_coord = 'none'
2291          mask_field%header%vertical_level = 1
2292
2293          call storage_get_field(mask_field, interp_mask_status)
2294
2295       end if
2296       if (interp_land_mask(idx) /= ' ') then
2297          mask_land_field%header%version = 1
2298          mask_land_field%header%forecast_hour = 0.
2299          mask_land_field%header%field = trim(interp_land_mask(idx))//'.mask'
2300          mask_land_field%header%vertical_coord = 'none'
2301          mask_land_field%header%vertical_level = 1
2302
2303          call storage_get_field(mask_land_field, interp_land_mask_status)
2304
2305       end if
2306       if (interp_water_mask(idx) /= ' ') then
2307          mask_water_field%header%version = 1
2308          mask_water_field%header%forecast_hour = 0.
2309          mask_water_field%header%field = trim(interp_water_mask(idx))//'.mask'
2310          mask_water_field%header%vertical_coord = 'none'
2311          mask_water_field%header%vertical_level = 1
2312
2313          call storage_get_field(mask_water_field, interp_water_mask_status)
2314
2315       end if
2316
2317       interp_array => interp_array_from_string(interp_method(idx))
2318       interp_opts => interp_options_from_string(interp_method(idx))
2319
2320
2321       !
2322       ! Interpolate using average_gcell interpolation method
2323       !
2324       if (do_gcell_interp) then
2325          !BPR BEGIN
2326          !If a lower priority source of the current field has already been read
2327          !in, the results of that input are currently in field%r_arr
2328          !Pass COPY of field%r_arr into accum_continous because in accum_continuous
2329          !it will set the input variable to zero over points covered by the
2330          !current source and then determine the appropriate value to place at
2331          !that point.  This will overwrite data already put in field%r_arr by
2332          !lower priority sources.
2333          !This is problematic if the current source results in missing values
2334          !over parts of the area covered by the current source where a lower
2335          !priority source has already provided a value. In this case, if one
2336          !passes in field%r_arr, it will overwrite the value provided by the
2337          !lower priority source with zero.
2338          r_arr_cur_source = field%r_arr
2339          !BPR END
2340          allocate(data_count(sm1:em1,sm2:em2))
2341          data_count = 0.
2342
2343          if (interp_mask_status == 0) then
2344             call accum_continuous(slab, &
2345                          minx, maxx, miny, maxy, 1, 1, bdr, &
2346                          !BPR BEGIN
2347                          !Pass COPY of field%r_arr instead of field%r_arr itself
2348                          !field%r_arr, data_count, &
2349                          r_arr_cur_source, data_count, &
2350                          !BPR END
2351                          sm1, em1, sm2, em2, 1, 1, &
2352                          istagger, &
2353                          new_pts, missing_value(idx), interp_mask_val(idx), interp_mask_relational(idx), mask_field%r_arr)
2354          else
2355             call accum_continuous(slab, &
2356                          minx, maxx, miny, maxy, 1, 1, bdr, &
2357                          !BPR BEGIN
2358                          !Pass COPY of field%r_arr instead of field%r_arr itself
2359                          !field%r_arr, data_count, &
2360                          r_arr_cur_source, data_count, &
2361                          !BPR END
2362                          sm1, em1, sm2, em2, 1, 1, &
2363                          istagger, &
2364                          new_pts, missing_value(idx), -1.) ! The -1 is the maskval, but since we
2365                                                            !   we do not give an optional mask, no
2366                                                            !   no need to worry about -1s in data
2367          end if
2368
2369          orig_selected_proj = iget_selected_domain()
2370          call select_domain(SOURCE_PROJ)
2371          do j=sm2,em2
2372             do i=sm1,em1
2373
2374                if (abs(i - sd1) >= process_width .and. (abs(i - ed1) >= process_width) .and. &
2375                    abs(j - sd2) >= process_width .and. (abs(j - ed2) >= process_width)) then
2376                   field%r_arr(i,j) = fill_missing(idx)
2377                   call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2378                   cycle
2379                end if
2380
2381                if (present(landmask)) then
2382
2383                   if (landmask(i,j) /= masked(idx)) then
2384                      if (data_count(i,j) > 0.) then
2385                         !BPR BEGIN
2386                         !accum_continuous is now passed a copy of field%r_arr (r_arr_cur_source)
2387                         !instead of field%r_arr itself so that it does not set
2388                         !field%r_arr to zero where the input source is marked as missing
2389                         !field%r_arr(i,j) = field%r_arr(i,j) / data_count(i,j)
2390                         field%r_arr(i,j) = r_arr_cur_source(i,j) / data_count(i,j)
2391                         !BPR END
2392                         call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2393                      else
2394
2395                         if (interp_mask_status == 0) then
2396                            temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2397                                                    minx, maxx, miny, maxy, bdr, missing_value(idx), &
2398                                                    mask_relational=interp_mask_relational(idx), &
2399                                                    mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
2400                         else
2401                            temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2402                                                    minx, maxx, miny, maxy, bdr, missing_value(idx))
2403                         end if
2404
2405                         if (temp /= missing_value(idx)) then
2406                            field%r_arr(i,j) = temp
2407                            call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2408                         end if
2409
2410                      end if
2411                   else
2412                      field%r_arr(i,j) = fill_missing(idx)
2413                      call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2414                   end if
2415
2416                   if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
2417                       .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
2418                      field%r_arr(i,j) = fill_missing(idx)
2419
2420                      ! Assume that if missing fill value is other than default, then user has asked
2421                      !    to fill in any missing values, and we can consider this point to have
2422                      !    received a valid value
2423                      if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2424                   end if
2425
2426                else
2427
2428                   if (data_count(i,j) > 0.) then
2429                      !BPR BEGIN
2430                      !accum_continuous is now passed a copy of field%r_arr (r_arr_cur_source)
2431                      !instead of field%r_arr itself so that it does not set
2432                      !field%r_arr to zero where the input source is marked as missing
2433                      !field%r_arr(i,j) = field%r_arr(i,j) / data_count(i,j)
2434                      field%r_arr(i,j) = r_arr_cur_source(i,j) / data_count(i,j)
2435                      !BPR END
2436                      call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2437                   else
2438
2439                      if (interp_mask_status == 0) then
2440                         temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2441                                                 minx, maxx, miny, maxy, bdr, missing_value(idx), &
2442                                                 mask_relational=interp_mask_relational(idx), &
2443                                                 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
2444                      else
2445                         temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2446                                                 minx, maxx, miny, maxy, bdr, missing_value(idx))
2447                      end if
2448
2449                      if (temp /= missing_value(idx)) then
2450                         field%r_arr(i,j) = temp
2451                         call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2452                      end if
2453
2454                   end if
2455
2456                   if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
2457                       .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
2458                      field%r_arr(i,j) = fill_missing(idx)
2459
2460                      ! Assume that if missing fill value is other than default, then user has asked
2461                      !    to fill in any missing values, and we can consider this point to have
2462                      !    received a valid value
2463                      if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2464                   end if
2465
2466                end if
2467
2468             end do
2469          end do
2470          call select_domain(orig_selected_proj)
2471          deallocate(data_count)
2472
2473       !
2474       ! No average_gcell interpolation method
2475       !
2476       else
2477
2478          orig_selected_proj = iget_selected_domain()
2479          call select_domain(SOURCE_PROJ)
2480          do j=sm2,em2
2481             do i=sm1,em1
2482
2483                if (abs(i - sd1) >= process_width .and. (abs(i - ed1) >= process_width) .and. &
2484                    abs(j - sd2) >= process_width .and. (abs(j - ed2) >= process_width)) then
2485                   field%r_arr(i,j) = fill_missing(idx)
2486                   call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2487                   cycle
2488                end if
2489
2490                if (present(landmask)) then
2491
2492                   if (masked(idx) == MASKED_BOTH) then
2493
2494                      if (landmask(i,j) == 0) then  ! WATER POINT
2495
2496                         if (interp_land_mask_status == 0) then
2497                            temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2498                                                    minx, maxx, miny, maxy, bdr, missing_value(idx), &
2499                                                    mask_relational=interp_land_mask_relational(idx), &
2500                                                    mask_val=interp_land_mask_val(idx), mask_field=mask_land_field%r_arr)
2501                         else
2502                            temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2503                                                    minx, maxx, miny, maxy, bdr, missing_value(idx))
2504                         end if
2505
2506                      else if (landmask(i,j) == 1) then  ! LAND POINT
2507
2508                         if (interp_water_mask_status == 0) then
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                                                    mask_relational=interp_water_mask_relational(idx), &
2512                                                    mask_val=interp_water_mask_val(idx), mask_field=mask_water_field%r_arr)
2513                         else
2514                            temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2515                                                    minx, maxx, miny, maxy, bdr, missing_value(idx))
2516                         end if
2517
2518                      end if
2519
2520                   else if (landmask(i,j) /= masked(idx)) then
2521
2522                      if (interp_mask_status == 0) then
2523                         temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2524                                                 minx, maxx, miny, maxy, bdr, missing_value(idx), &
2525                                                 mask_relational=interp_mask_relational(idx), &
2526                                                 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
2527                      else
2528                         temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2529                                                 minx, maxx, miny, maxy, bdr, missing_value(idx))
2530                      end if
2531
2532                   else
2533                      temp = missing_value(idx)
2534                   end if
2535
2536                ! No landmask for this field
2537                else
2538
2539                   if (interp_mask_status == 0) then
2540                      temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2541                                              minx, maxx, miny, maxy, bdr, missing_value(idx), &
2542                                              mask_relational=interp_mask_relational(idx), &
2543                                              mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
2544                   else
2545                      temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
2546                                              minx, maxx, miny, maxy, bdr, missing_value(idx))
2547                   end if
2548
2549                end if
2550
2551                if (temp /= missing_value(idx)) then
2552                   field%r_arr(i,j) = temp
2553                   call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2554                else if (present(landmask)) then
2555                   if (landmask(i,j) == masked(idx)) then
2556                      field%r_arr(i,j) = fill_missing(idx)
2557                      call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2558                   end if
2559                end if
2560
2561                if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
2562                    .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
2563                   field%r_arr(i,j) = fill_missing(idx)
2564
2565                   ! Assume that if missing fill value is other than default, then user has asked
2566                   !    to fill in any missing values, and we can consider this point to have
2567                   !    received a valid value
2568                   if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
2569                end if
2570
2571             end do
2572          end do
2573          call select_domain(orig_selected_proj)
2574       end if
2575
2576       deallocate(interp_array)
2577       deallocate(interp_opts)
2578
2579    end subroutine interp_met_field
2580
2581
2582    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2583    ! Name: interp_to_latlon
2584    !
2585    ! Purpose:
2586    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2587    function interp_to_latlon(rlat, rlon, istagger, interp_method_list, interp_opt_list, slab, &
2588                              minx, maxx, miny, maxy, bdr, source_missing_value, &
2589                              mask_field, mask_relational, mask_val)
2590
2591       use interp_module
2592       use llxy_module
2593
2594       implicit none
2595
2596       ! Arguments
2597       integer, intent(in) :: minx, maxx, miny, maxy, bdr, istagger
2598       integer, dimension(:), intent(in) :: interp_method_list
2599       integer, dimension(:), intent(in) :: interp_opt_list
2600       real, intent(in) :: rlat, rlon, source_missing_value
2601       real, dimension(minx:maxx,miny:maxy), intent(in) :: slab
2602       real, intent(in), optional :: mask_val
2603       real, dimension(minx:maxx,miny:maxy), intent(in), optional :: mask_field
2604       character(len=1), intent(in), optional :: mask_relational
2605
2606       ! Return value
2607       real :: interp_to_latlon
2608
2609       ! Local variables
2610       real :: rx, ry
2611
2612       interp_to_latlon = source_missing_value
2613
2614       call lltoxy(rlat, rlon, rx, ry, istagger)
2615       if (rx >= minx+bdr-0.5 .and. rx <= maxx-bdr+0.5) then
2616          if (present(mask_field) .and. present(mask_val) .and. present (mask_relational)) then
2617             interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
2618                                                interp_method_list, interp_opt_list, 1, mask_relational, mask_val, mask_field)
2619          else if (present(mask_field) .and. present(mask_val)) then
2620             interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
2621                                                interp_method_list, interp_opt_list, 1, maskval=mask_val, mask_array=mask_field)
2622          else
2623             interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
2624                                                interp_method_list, interp_opt_list, 1)
2625          end if
2626       else
2627          interp_to_latlon = source_missing_value
2628       end if
2629
2630       if (interp_to_latlon == source_missing_value) then
2631
2632          ! Try a lon in the range 0. to 360.; all lons in the xlon
2633          !    array should be in the range -180. to 180.
2634          if (rlon < 0.) then
2635             call lltoxy(rlat, rlon+360., rx, ry, istagger)
2636             if (rx >= minx+bdr-0.5 .and. rx <= maxx-bdr+0.5) then
2637                if (present(mask_field) .and. present(mask_val) .and. present(mask_relational)) then
2638                   interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
2639                                                      1, 1, source_missing_value, &
2640                                                      interp_method_list, interp_opt_list, 1, &
2641                                                      mask_relational, mask_val, mask_field)
2642                else if (present(mask_field) .and. present(mask_val)) then
2643                   interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
2644                                                      1, 1, source_missing_value, &
2645                                                      interp_method_list, interp_opt_list, 1, &
2646                                                      maskval=mask_val, mask_array=mask_field)
2647                else
2648                   interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
2649                                                      1, 1, source_missing_value, &
2650                                                      interp_method_list, interp_opt_list, 1)
2651                end if
2652             else
2653                interp_to_latlon = source_missing_value
2654             end if
2655
2656          end if
2657
2658       end if
2659
2660       return
2661
2662    end function interp_to_latlon
2663
2664
2665    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2666    ! Name: get_bottom_top_dim
2667    !
2668    ! Purpose:
2669    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2670    subroutine get_bottom_top_dim(bottom_top_dim)
2671
2672       use interp_option_module
2673       use list_module
2674       use storage_module
2675
2676       implicit none
2677
2678       ! Arguments
2679       integer, intent(out) :: bottom_top_dim
2680
2681       ! Local variables
2682       integer :: i, j
2683       integer, pointer, dimension(:) :: field_levels
2684       character (len=32) :: z_dim
2685       type (fg_input), pointer, dimension(:) :: headers
2686       type (list) :: temp_levels
2687
2688       !CWH Initialize local pointer variables
2689       nullify(field_levels)
2690       nullify(headers)
2691
2692       ! Initialize a list to store levels that are found for 3-d fields
2693       call list_init(temp_levels)
2694
2695       ! Get a list of all time-dependent fields (given by their headers) from
2696       !   the storage module
2697       call storage_get_td_headers(headers)
2698
2699       !
2700       ! Given headers of all fields, we first build a list of all possible levels
2701       !    for 3-d met fields (excluding sea-level, though).
2702       !
2703       do i=1,size(headers)
2704          call get_z_dim_name(headers(i)%header%field, z_dim)
2705
2706          ! We only want to consider 3-d met fields
2707          if (z_dim(1:18) == 'num_metgrid_levels') then
2708
2709             ! Find out what levels the current field has
2710             call storage_get_levels(headers(i), field_levels)
2711             do j=1,size(field_levels)
2712
2713                ! If this level has not yet been encountered, add it to our list
2714                if (.not. list_search(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))) then
2715                   if (field_levels(j) /= 201300) then
2716                      call list_insert(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))
2717                   end if
2718                end if
2719
2720             end do
2721
2722             deallocate(field_levels)
2723
2724          end if
2725
2726       end do
2727
2728       bottom_top_dim = list_length(temp_levels)
2729
2730       call list_destroy(temp_levels)
2731       deallocate(headers)
2732
2733    end subroutine get_bottom_top_dim
2734
2735
2736    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2737    ! Name: fill_missing_levels
2738    !
2739    ! Purpose:
2740    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2741    subroutine fill_missing_levels(output_flags)
2742
2743       use interp_option_module
2744       use list_module
2745       use module_debug
2746       use module_mergesort
2747       use storage_module
2748
2749       implicit none
2750
2751       ! Arguments
2752       character (len=128), dimension(:), intent(inout) :: output_flags
2753
2754       ! Local variables
2755       integer :: i, ii, j, ix, jx, k, lower, upper, temp, istatus
2756       integer, pointer, dimension(:) :: union_levels, field_levels
2757       real, pointer, dimension(:) :: r_union_levels
2758       character (len=128) :: clevel
2759       type (fg_input) :: lower_field, upper_field, new_field, search_field
2760       type (fg_input), pointer, dimension(:) :: headers, all_headers
2761       type (list) :: temp_levels
2762       type (list_item), pointer, dimension(:) :: keys
2763
2764       ! CWH Initialize local pointer variables
2765       nullify(union_levels)
2766       nullify(field_levels)
2767       nullify(r_union_levels)
2768       nullify(headers)
2769       nullify(all_headers)
2770       nullify(keys)
2771
2772       ! Initialize a list to store levels that are found for 3-d fields
2773       call list_init(temp_levels)
2774
2775       ! Get a list of all fields (given by their headers) from the storage module
2776       call storage_get_td_headers(headers)
2777       call storage_get_all_headers(all_headers)
2778
2779       !
2780       ! Given headers of all fields, we first build a list of all possible levels
2781       !    for 3-d met fields (excluding sea-level, though).
2782       !
2783       do i=1,size(headers)
2784
2785          ! Find out what levels the current field has
2786          call storage_get_levels(headers(i), field_levels)
2787          do j=1,size(field_levels)
2788
2789             ! If this level has not yet been encountered, add it to our list
2790             if (.not. list_search(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))) then
2791                if (field_levels(j) /= 201300) then
2792                   call list_insert(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))
2793                end if
2794             end if
2795
2796          end do
2797
2798          deallocate(field_levels)
2799
2800       end do
2801
2802       if (list_length(temp_levels) > 0) then
2803
2804          !
2805          ! With all possible levels stored in a list, get an array of levels, sorted
2806          !    in decreasing order
2807          !
2808          i = 0
2809          allocate(union_levels(list_length(temp_levels)))
2810          do while (list_length(temp_levels) > 0)
2811             i = i + 1
2812             call list_get_first_item(temp_levels, ikey=union_levels(i), ivalue=temp)
2813          end do
2814          call mergesort(union_levels, 1, size(union_levels))
2815
2816          allocate(r_union_levels(size(union_levels)))
2817          do i=1,size(union_levels)
2818             r_union_levels(i) = real(union_levels(i))
2819          end do
2820
2821          !
2822          ! With a sorted, complete list of levels, we need
2823          !    to go back and fill in missing levels for each 3-d field
2824          !
2825          do i=1,size(headers)
2826
2827             !
2828             ! Find entry in METGRID.TBL for this field, if one exists; if it does, then the
2829             !    entry may tell us how to get values for the current field at the missing level
2830             !
2831             do ii=1,num_entries
2832                if (fieldname(ii) == headers(i)%header%field) exit
2833             end do
2834             if (ii <= num_entries) then
2835                call dup(headers(i),new_field)
2836                nullify(new_field%valid_mask)
2837                nullify(new_field%modified_mask)
2838                call fill_field(new_field, ii, output_flags, r_union_levels)
2839             end if
2840
2841          end do
2842
2843          deallocate(union_levels)
2844          deallocate(r_union_levels)
2845          deallocate(headers)
2846
2847          call storage_get_td_headers(headers)
2848
2849          !
2850          ! Now we may need to vertically interpolate to missing values in 3-d fields
2851          !
2852          do i=1,size(headers)
2853
2854             call storage_get_levels(headers(i), field_levels)
2855
2856             ! If this isn't a 3-d array, nothing to do
2857             if (size(field_levels) > 1) then
2858
2859                do k=1,size(field_levels)
2860                   call dup(headers(i),search_field)
2861                   search_field%header%vertical_level = field_levels(k)
2862                   call storage_get_field(search_field,istatus)
2863                   if (istatus == 0) then
2864                      JLOOP: do jx=search_field%header%dim2(1),search_field%header%dim2(2)
2865                         ILOOP: do ix=search_field%header%dim1(1),search_field%header%dim1(2)
2866                            if (.not. bitarray_test(search_field%valid_mask, &
2867                                                    ix-search_field%header%dim1(1)+1, &
2868                                                    jx-search_field%header%dim2(1)+1)) then
2869
2870                               call dup(search_field, lower_field)
2871                               do lower=k-1,1,-1
2872                                  lower_field%header%vertical_level = field_levels(lower)
2873                                  call storage_get_field(lower_field,istatus)
2874                                  if (bitarray_test(lower_field%valid_mask, &
2875                                                    ix-search_field%header%dim1(1)+1, &
2876                                                    jx-search_field%header%dim2(1)+1)) &
2877                                      exit
2878
2879                               end do
2880
2881                               call dup(search_field, upper_field)
2882                               do upper=k+1,size(field_levels)
2883                                  upper_field%header%vertical_level = field_levels(upper)
2884                                  call storage_get_field(upper_field,istatus)
2885                                  if (bitarray_test(upper_field%valid_mask, &
2886                                                    ix-search_field%header%dim1(1)+1, &
2887                                                    jx-search_field%header%dim2(1)+1)) &
2888                                      exit
2889
2890                               end do
2891                               if (upper <= size(field_levels) .and. lower >= 1) then
2892                                  search_field%r_arr(ix,jx) = real(abs(field_levels(upper)-field_levels(k))) &
2893                                                            / real(abs(field_levels(upper)-field_levels(lower))) &
2894                                                            * lower_field%r_arr(ix,jx) &
2895                                                            + real(abs(field_levels(k)-field_levels(lower))) &
2896                                                            / real(abs(field_levels(upper)-field_levels(lower))) &
2897                                                            * upper_field%r_arr(ix,jx)
2898                                  call bitarray_set(search_field%valid_mask, &
2899                                                    ix-search_field%header%dim1(1)+1, &
2900                                                    jx-search_field%header%dim2(1)+1)
2901                               end if
2902                            end if
2903                         end do ILOOP
2904                      end do JLOOP
2905                   else
2906                      call mprintf(.true.,ERROR, &
2907                                   'This is bad, could not get %s at level %i.', &
2908                                   s1=trim(search_field%header%field), i1=field_levels(k))
2909                   end if
2910                end do
2911
2912             end if
2913
2914             deallocate(field_levels)
2915
2916          end do
2917
2918       end if
2919
2920       call list_destroy(temp_levels)
2921       deallocate(all_headers)
2922       deallocate(headers)
2923
2924    end subroutine fill_missing_levels
2925
2926
2927    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2928    ! Name: create_derived_fields
2929    !
2930    ! Purpose:
2931    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2932    subroutine create_derived_fields(arg_gridtype, hdate, xfcst, &
2933                                  we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
2934                                  we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e, &
2935                                  created_this_field, output_flags)
2936
2937       use interp_option_module
2938       use list_module
2939       use module_mergesort
2940       use storage_module
2941
2942       implicit none
2943
2944       ! Arguments
2945       integer, intent(in) :: we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
2946                              we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e
2947       real, intent(in) :: xfcst
2948       logical, dimension(:), intent(inout) :: created_this_field
2949       character (len=1), intent(in) :: arg_gridtype
2950       character (len=24), intent(in) :: hdate
2951       character (len=128), dimension(:), intent(inout) :: output_flags
2952
2953       ! Local variables
2954       integer :: idx, i, j, istatus
2955       type (fg_input) :: field
2956
2957       ! Initialize fg_input structure to store the field
2958       field%header%version = 1
2959       field%header%date = hdate//'        '
2960       field%header%time_dependent = .true.
2961       field%header%mask_field = .false.
2962       field%header%constant_field = .false.
2963       field%header%forecast_hour = xfcst
2964       field%header%fg_source = 'Derived from FG'
2965       field%header%field = ' '
2966       field%header%units = ' '
2967       field%header%description = ' '
2968       field%header%vertical_level = 0
2969       field%header%sr_x = 1
2970       field%header%sr_y = 1
2971       field%header%array_order = 'XY '
2972       field%header%is_wind_grid_rel = .true.
2973       field%header%array_has_missing_values = .false.
2974       nullify(field%r_arr)
2975       nullify(field%valid_mask)
2976       nullify(field%modified_mask)
2977
2978       !
2979       ! Check each entry in METGRID.TBL to see whether it is a derive field
2980       !
2981       do idx=1,num_entries
2982          if (is_derived_field(idx)) then
2983
2984             created_this_field(idx) = .true.
2985
2986             call mprintf(.true.,INFORM,'Going to create the field %s',s1=fieldname(idx))
2987
2988             ! Intialize more fields in storage structure
2989             field%header%field = fieldname(idx)
2990             call get_z_dim_name(fieldname(idx),field%header%vertical_coord)
2991             field%map%stagger = output_stagger(idx)
2992             if (arg_gridtype == 'E') then
2993                field%header%dim1(1) = we_mem_s
2994                field%header%dim1(2) = we_mem_e
2995                field%header%dim2(1) = sn_mem_s
2996                field%header%dim2(2) = sn_mem_e
2997             else if (arg_gridtype == 'C') then
2998                if (output_stagger(idx) == M) then
2999                   field%header%dim1(1) = we_mem_s
3000                   field%header%dim1(2) = we_mem_e
3001                   field%header%dim2(1) = sn_mem_s
3002                   field%header%dim2(2) = sn_mem_e
3003                else if (output_stagger(idx) == U) then
3004                   field%header%dim1(1) = we_mem_stag_s
3005                   field%header%dim1(2) = we_mem_stag_e
3006                   field%header%dim2(1) = sn_mem_s
3007                   field%header%dim2(2) = sn_mem_e
3008                else if (output_stagger(idx) == V) then
3009                   field%header%dim1(1) = we_mem_s
3010                   field%header%dim1(2) = we_mem_e
3011                   field%header%dim2(1) = sn_mem_stag_s
3012                   field%header%dim2(2) = sn_mem_stag_e
3013                end if
3014             end if
3015
3016             call fill_field(field, idx, output_flags)
3017
3018          end if
3019       end do
3020
3021
3022    end subroutine create_derived_fields
3023
3024
3025    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3026    ! Name: fill_field
3027    !
3028    ! Purpose:
3029    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3030    subroutine fill_field(field, idx, output_flags, all_level_list)
3031
3032       use interp_option_module
3033       use list_module
3034       use module_mergesort
3035       use storage_module
3036
3037       implicit none
3038
3039       ! Arguments
3040       integer, intent(in) :: idx
3041       type (fg_input), intent(inout) :: field
3042       character (len=128), dimension(:), intent(inout) :: output_flags
3043       real, dimension(:), intent(in), optional :: all_level_list
3044
3045       ! Local variables
3046       integer :: i, j, istatus, isrclevel
3047       integer, pointer, dimension(:) :: all_list
3048       real :: rfillconst, rlevel, rsrclevel
3049       type (fg_input) :: query_field
3050       type (list_item), pointer, dimension(:) :: keys
3051       character (len=128) :: asrcname
3052       logical :: filled_all_lev
3053
3054       !CWH Initialize local pointer variables
3055       nullify(all_list)
3056       nullify(keys)
3057
3058       filled_all_lev = .false.
3059
3060       !
3061       ! Get a list of all levels to be filled for this field
3062       !
3063       keys => list_get_keys(fill_lev_list(idx))
3064
3065       do i=1,list_length(fill_lev_list(idx))
3066
3067          !
3068          ! First handle a specification for levels "all"
3069          !
3070          if (trim(keys(i)%ckey) == 'all') then
3071
3072             ! We only want to fill all levels if we haven't already filled "all" of them
3073             if (.not. filled_all_lev) then
3074
3075                filled_all_lev = .true.
3076
3077                query_field%header%time_dependent = .true.
3078                query_field%header%field = ' '
3079                nullify(query_field%r_arr)
3080                nullify(query_field%valid_mask)
3081                nullify(query_field%modified_mask)
3082
3083                ! See if we are filling this level with a constant
3084                call get_constant_fill_lev(keys(i)%cvalue, rfillconst, istatus)
3085                if (istatus == 0) then
3086                   if (present(all_level_list)) then
3087                      do j=1,size(all_level_list)
3088                         call create_level(field, real(all_level_list(j)), idx, output_flags, rfillconst=rfillconst)
3089                      end do
3090                   else
3091                      query_field%header%field = level_template(idx)
3092                      nullify(all_list)
3093                      call storage_get_levels(query_field, all_list)
3094                      if (associated(all_list)) then
3095                         do j=1,size(all_list)
3096                            call create_level(field, real(all_list(j)), idx, output_flags, rfillconst=rfillconst)
3097                         end do
3098                         deallocate(all_list)
3099                      end if
3100                   end if
3101
3102                ! Else see if we are filling this level with a constant equal
3103                !   to the value of the level
3104                else if (trim(keys(i)%cvalue) == 'vertical_index') then
3105                   if (present(all_level_list)) then
3106                      do j=1,size(all_level_list)
3107                         call create_level(field, real(all_level_list(j)), idx, output_flags, &
3108                                           rfillconst=real(all_level_list(j)))
3109                      end do
3110                   else
3111                      query_field%header%field = level_template(idx)
3112                      nullify(all_list)
3113                      call storage_get_levels(query_field, all_list)
3114                      if (associated(all_list)) then
3115                         do j=1,size(all_list)
3116                            call create_level(field, real(all_list(j)), idx, output_flags, rfillconst=real(all_list(j)))
3117                         end do
3118                         deallocate(all_list)
3119                      end if
3120                   end if
3121
3122                ! Else, we assume that it is a field from which we are copying levels
3123                else
3124                   if (present(all_level_list)) then
3125                      do j=1,size(all_level_list)
3126                         call create_level(field, real(all_level_list(j)), idx, output_flags, &
3127                                           asrcname=keys(i)%cvalue, rsrclevel=real(all_level_list(j)))
3128                      end do
3129                   else
3130                      query_field%header%field = keys(i)%cvalue  ! Use same levels as source field, not level_template
3131                      nullify(all_list)
3132                      call storage_get_levels(query_field, all_list)
3133                      if (associated(all_list)) then
3134                         do j=1,size(all_list)
3135                            call create_level(field, real(all_list(j)), idx, output_flags, &
3136                                              asrcname=keys(i)%cvalue, rsrclevel=real(all_list(j)))
3137                         end do
3138                         deallocate(all_list)
3139
3140                      else
3141
3142                         ! If the field doesn't have any levels (or does not exist) then we have not
3143                         !   really filled all levels at this point.
3144                         filled_all_lev = .false.
3145                      end if
3146                   end if
3147
3148                end if
3149             end if
3150
3151          !
3152          ! Handle individually specified levels
3153          !
3154          else
3155
3156             read(keys(i)%ckey,*) rlevel
3157
3158             ! See if we are filling this level with a constant
3159             call get_constant_fill_lev(keys(i)%cvalue, rfillconst, istatus)
3160             if (istatus == 0) then
3161                call create_level(field, rlevel, idx, output_flags, rfillconst=rfillconst)
3162
3163             ! Otherwise, we are filling from another level
3164             else
3165                call get_fill_src_level(keys(i)%cvalue, asrcname, isrclevel)
3166                rsrclevel = real(isrclevel)
3167                call create_level(field, rlevel, idx, output_flags, &
3168                                  asrcname=asrcname, rsrclevel=rsrclevel)
3169
3170             end if
3171          end if
3172       end do
3173
3174       if (associated(keys)) deallocate(keys)
3175
3176    end subroutine fill_field
3177
3178
3179    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3180    ! Name: create_level
3181    !
3182    ! Purpose:
3183    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3184    subroutine create_level(field_template, rlevel, idx, output_flags, &
3185                            rfillconst, asrcname, rsrclevel)
3186
3187       use storage_module
3188       use interp_option_module
3189
3190       implicit none
3191
3192       ! Arguments
3193       type (fg_input), intent(inout) :: field_template
3194       real, intent(in) :: rlevel
3195       integer, intent(in) :: idx
3196       character (len=128), dimension(:), intent(inout) :: output_flags
3197       real, intent(in), optional :: rfillconst, rsrclevel
3198       character (len=128), intent(in), optional :: asrcname
3199
3200       ! Local variables
3201       integer :: i, j, istatus
3202       integer :: sm1,em1,sm2,em2
3203       type (fg_input) :: query_field
3204
3205       !
3206       ! Check to make sure optional arguments are sane
3207       !
3208       if (present(rfillconst) .and. (present(asrcname) .or. present(rsrclevel))) then
3209          call mprintf(.true.,ERROR,'A call to create_level() cannot be given specifications '// &
3210                       'for both a constant fill value and a source level.')
3211
3212       else if ((present(asrcname) .and. .not. present(rsrclevel)) .or. &
3213                (.not. present(asrcname) .and. present(rsrclevel))) then
3214          call mprintf(.true.,ERROR,'Neither or both of optional arguments asrcname and '// &
3215                       'rsrclevel must be specified to subroutine create_level().')
3216
3217       else if (.not. present(rfillconst) .and. &
3218                .not. present(asrcname)   .and. &
3219                .not. present(rsrclevel)) then
3220          call mprintf(.true.,ERROR,'A call to create_level() must be given either a specification '// &
3221                       'for a constant fill value or a source level.')
3222       end if
3223
3224       query_field%header%time_dependent = .true.
3225       query_field%header%field = field_template%header%field
3226       query_field%header%vertical_level = rlevel
3227       nullify(query_field%r_arr)
3228       nullify(query_field%valid_mask)
3229       nullify(query_field%modified_mask)
3230
3231       call storage_query_field(query_field, istatus)
3232       if (istatus == 0) then
3233          call mprintf(.true.,INFORM,'%s at level %f already exists; leaving it alone.', &
3234                       s1=field_template%header%field, f1=rlevel)
3235          return
3236       end if
3237
3238       sm1 = field_template%header%dim1(1)
3239       em1 = field_template%header%dim1(2)
3240       sm2 = field_template%header%dim2(1)
3241       em2 = field_template%header%dim2(2)
3242
3243       !
3244       ! Handle constant fill value case
3245       !
3246       if (present(rfillconst)) then
3247
3248          field_template%header%vertical_level = rlevel
3249          allocate(field_template%r_arr(sm1:em1,sm2:em2))
3250          allocate(field_template%valid_mask)
3251          allocate(field_template%modified_mask)
3252          call bitarray_create(field_template%valid_mask, em1-sm1+1, em2-sm2+1)
3253          call bitarray_create(field_template%modified_mask, em1-sm1+1, em2-sm2+1)
3254
3255          field_template%r_arr = rfillconst
3256
3257          do j=sm2,em2
3258             do i=sm1,em1
3259                call bitarray_set(field_template%valid_mask, i-sm1+1, j-sm2+1)
3260             end do
3261          end do
3262
3263          call storage_put_field(field_template)
3264
3265          if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
3266             output_flags(idx) = flag_in_output(idx)
3267          end if
3268
3269       !
3270       ! Handle source field and source level case
3271       !
3272       else if (present(asrcname) .and. present(rsrclevel)) then
3273
3274          query_field%header%field = ' '
3275          query_field%header%field = asrcname
3276          query_field%header%vertical_level = rsrclevel
3277
3278          ! Check to see whether the requested source field exists at the requested level
3279          call storage_query_field(query_field, istatus)
3280
3281          if (istatus == 0) then
3282
3283             ! Read in requested field at requested level
3284             call storage_get_field(query_field, istatus)
3285             if ((query_field%header%dim1(1) /= field_template%header%dim1(1)) .or. &
3286                 (query_field%header%dim1(2) /= field_template%header%dim1(2)) .or. &
3287                 (query_field%header%dim2(1) /= field_template%header%dim2(1)) .or. &
3288                 (query_field%header%dim2(2) /= field_template%header%dim2(2))) then
3289                call mprintf(.true.,ERROR,'Dimensions for %s do not match those of %s. This is '// &
3290                             'probably because the staggerings of the fields do not match.', &
3291                             s1=query_field%header%field, s2=field_template%header%field)
3292             end if
3293
3294             field_template%header%vertical_level = rlevel
3295             allocate(field_template%r_arr(sm1:em1,sm2:em2))
3296             allocate(field_template%valid_mask)
3297             allocate(field_template%modified_mask)
3298             call bitarray_create(field_template%valid_mask, em1-sm1+1, em2-sm2+1)
3299             call bitarray_create(field_template%modified_mask, em1-sm1+1, em2-sm2+1)
3300
3301             field_template%r_arr = query_field%r_arr
3302
3303             ! We should retain information about which points in the field are valid
3304             do j=sm2,em2
3305                do i=sm1,em1
3306                   if (bitarray_test(query_field%valid_mask, i-sm1+1, j-sm2+1)) then
3307                      call bitarray_set(field_template%valid_mask, i-sm1+1, j-sm2+1)
3308                   end if
3309                end do
3310             end do
3311
3312             call storage_put_field(field_template)
3313
3314             if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
3315                output_flags(idx) = flag_in_output(idx)
3316             end if
3317
3318          else
3319             call mprintf(.true.,INFORM,'Couldn''t find %s at level %f to fill level %f of %s.', &
3320                          s1=asrcname,f1=rsrclevel,f2=rlevel,s2=field_template%header%field)
3321          end if
3322
3323       end if
3324
3325    end subroutine create_level
3326
3327
3328    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3329    ! Name: accum_continuous
3330    !
3331    ! Purpose: Sum up all of the source data points whose nearest neighbor in the
3332    !   model grid is the specified model grid point.
3333    !
3334    ! NOTE: When processing the source tile, those source points that are
3335    !   closest to a different model grid point will be added to the totals for
3336    !   such grid points; thus, an entire source tile will be processed at a time.
3337    !   This routine really processes for all model grid points that are
3338    !   within a source tile, and not just for a single grid point.
3339    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3340    subroutine accum_continuous(src_array, &
3341                                src_min_x, src_max_x, src_min_y, src_max_y, src_min_z, src_max_z, bdr_width, &
3342                                dst_array, n, &
3343                                start_i, end_i, start_j, end_j, start_k, end_k, &
3344                                istagger, &
3345                                new_pts, msgval, maskval, mask_relational, mask_array, sr_x, sr_y)
3346
3347       use bitarray_module
3348       use misc_definitions_module
3349
3350       implicit none
3351
3352       ! Arguments
3353       integer, intent(in) :: start_i, end_i, start_j, end_j, start_k, end_k, istagger, &
3354                              src_min_x, src_max_x, src_min_y, src_max_y, src_min_z, src_max_z, bdr_width
3355       real, intent(in) :: maskval, msgval
3356       real, dimension(src_min_x:src_max_x, src_min_y:src_max_y, src_min_z:src_max_z), intent(in) :: src_array
3357       real, dimension(start_i:end_i, start_j:end_j, start_k:end_k), intent(inout) :: dst_array, n
3358       real, dimension(src_min_x:src_max_x, src_min_y:src_max_y), intent(in), optional :: mask_array
3359       integer, intent(in), optional :: sr_x, sr_y
3360       type (bitarray), intent(inout) :: new_pts
3361       character(len=1), intent(in), optional :: mask_relational
3362
3363       ! Local variables
3364       integer :: i, j
3365       integer, pointer, dimension(:,:,:) :: where_maps_to
3366       real :: rsr_x, rsr_y
3367
3368       rsr_x = 1.0
3369       rsr_y = 1.0
3370       if (present(sr_x)) rsr_x = real(sr_x)
3371       if (present(sr_y)) rsr_y = real(sr_y)
3372
3373       allocate(where_maps_to(src_min_x:src_max_x,src_min_y:src_max_y,2))
3374       do i=src_min_x,src_max_x
3375          do j=src_min_y,src_max_y
3376             where_maps_to(i,j,1) = NOT_PROCESSED
3377          end do
3378       end do
3379
3380       call process_continuous_block(src_array, where_maps_to, &
3381                                src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
3382                                src_min_x+bdr_width, src_min_y, src_min_z, &
3383                                src_max_x-bdr_width, src_max_y, src_max_z, &
3384                                dst_array, n, start_i, end_i, start_j, end_j, start_k, end_k, &
3385                                istagger, &
3386                                new_pts, rsr_x, rsr_y, msgval, maskval, mask_relational, mask_array)
3387
3388       deallocate(where_maps_to)
3389
3390    end subroutine accum_continuous
3391
3392
3393    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3394    ! Name: process_continuous_block
3395    !
3396    ! Purpose: To recursively process a subarray of continuous data, adding the
3397    !   points in a block to the sum for their nearest grid point. The nearest
3398    !   neighbor may be estimated in some cases; for example, if the four corners
3399    !   of a subarray all have the same nearest grid point, all elements in the
3400    !   subarray are added to that grid point.
3401    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3402    recursive subroutine process_continuous_block(tile_array, where_maps_to, &
3403                                    src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
3404                                    min_i, min_j, min_k, max_i, max_j, max_k, &
3405                                    dst_array, n, &
3406                                    start_x, end_x, start_y, end_y, start_z, end_z, &
3407                                    istagger, &
3408                                    new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
3409
3410       use bitarray_module
3411       use llxy_module
3412       use misc_definitions_module
3413
3414       implicit none
3415
3416       ! Arguments
3417       integer, intent(in) :: min_i, min_j, min_k, max_i, max_j, max_k, &
3418                              src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
3419                              start_x, end_x, start_y, end_y, start_z, end_z, istagger
3420       integer, dimension(src_min_x:src_max_x,src_min_y:src_max_y,2), intent(inout) :: where_maps_to
3421       real, intent(in) :: sr_x, sr_y, maskval, msgval
3422       real, dimension(src_min_x:src_max_x,src_min_y:src_max_y,src_min_z:src_max_z), intent(in) :: tile_array
3423       real, dimension(src_min_x:src_max_x,src_min_y:src_max_y), intent(in), optional :: mask_array
3424       real, dimension(start_x:end_x,start_y:end_y,start_z:end_z), intent(inout) :: dst_array, n
3425       type (bitarray), intent(inout) :: new_pts
3426       character(len=1), intent(in), optional :: mask_relational
3427
3428       ! Local variables
3429       integer :: orig_selected_domain, x_dest, y_dest, i, j, k, center_i, center_j
3430       real :: lat_corner, lon_corner, rx, ry
3431
3432       ! Compute the model grid point that the corners of the rectangle to be
3433       !   processed map to
3434       ! Lower-left corner
3435       if (where_maps_to(min_i,min_j,1) == NOT_PROCESSED) then
3436          orig_selected_domain = iget_selected_domain()
3437          call select_domain(SOURCE_PROJ)
3438          call xytoll(real(min_i), real(min_j), lat_corner, lon_corner, istagger)
3439          call select_domain(1)
3440          call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
3441          rx = (rx - 1.0)*sr_x + 1.0
3442          ry = (ry - 1.0)*sr_y + 1.0
3443          call select_domain(orig_selected_domain)
3444          if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
3445              real(start_y) <= ry .and. ry <= real(end_y)) then
3446             where_maps_to(min_i,min_j,1) = nint(rx)
3447             where_maps_to(min_i,min_j,2) = nint(ry)
3448          else
3449             where_maps_to(min_i,min_j,1) = OUTSIDE_DOMAIN
3450          end if
3451       end if
3452
3453       ! Upper-left corner
3454       if (where_maps_to(min_i,max_j,1) == NOT_PROCESSED) then
3455          orig_selected_domain = iget_selected_domain()
3456          call select_domain(SOURCE_PROJ)
3457          call xytoll(real(min_i), real(max_j), lat_corner, lon_corner, istagger)
3458          call select_domain(1)
3459          call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
3460          rx = (rx - 1.0)*sr_x + 1.0
3461          ry = (ry - 1.0)*sr_y + 1.0
3462          call select_domain(orig_selected_domain)
3463          if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
3464              real(start_y) <= ry .and. ry <= real(end_y)) then
3465             where_maps_to(min_i,max_j,1) = nint(rx)
3466             where_maps_to(min_i,max_j,2) = nint(ry)
3467          else
3468             where_maps_to(min_i,max_j,1) = OUTSIDE_DOMAIN
3469          end if
3470       end if
3471
3472       ! Upper-right corner
3473       if (where_maps_to(max_i,max_j,1) == NOT_PROCESSED) then
3474          orig_selected_domain = iget_selected_domain()
3475          call select_domain(SOURCE_PROJ)
3476          call xytoll(real(max_i), real(max_j), lat_corner, lon_corner, istagger)
3477          call select_domain(1)
3478          call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
3479          rx = (rx - 1.0)*sr_x + 1.0
3480          ry = (ry - 1.0)*sr_y + 1.0
3481          call select_domain(orig_selected_domain)
3482          if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
3483              real(start_y) <= ry .and. ry <= real(end_y)) then
3484             where_maps_to(max_i,max_j,1) = nint(rx)
3485             where_maps_to(max_i,max_j,2) = nint(ry)
3486          else
3487             where_maps_to(max_i,max_j,1) = OUTSIDE_DOMAIN
3488          end if
3489       end if
3490
3491       ! Lower-right corner
3492       if (where_maps_to(max_i,min_j,1) == NOT_PROCESSED) then
3493          orig_selected_domain = iget_selected_domain()
3494          call select_domain(SOURCE_PROJ)
3495          call xytoll(real(max_i), real(min_j), lat_corner, lon_corner, istagger)
3496          call select_domain(1)
3497          call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
3498          rx = (rx - 1.0)*sr_x + 1.0
3499          ry = (ry - 1.0)*sr_y + 1.0
3500          call select_domain(orig_selected_domain)
3501          if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
3502              real(start_y) <= ry .and. ry <= real(end_y)) then
3503             where_maps_to(max_i,min_j,1) = nint(rx)
3504             where_maps_to(max_i,min_j,2) = nint(ry)
3505          else
3506             where_maps_to(max_i,min_j,1) = OUTSIDE_DOMAIN
3507          end if
3508       end if
3509
3510       ! If all four corners map to same model grid point, accumulate the
3511       !   entire rectangle
3512       if (where_maps_to(min_i,min_j,1) == where_maps_to(min_i,max_j,1) .and. &
3513           where_maps_to(min_i,min_j,1) == where_maps_to(max_i,max_j,1) .and. &
3514           where_maps_to(min_i,min_j,1) == where_maps_to(max_i,min_j,1) .and. &
3515           where_maps_to(min_i,min_j,2) == where_maps_to(min_i,max_j,2) .and. &
3516           where_maps_to(min_i,min_j,2) == where_maps_to(max_i,max_j,2) .and. &
3517           where_maps_to(min_i,min_j,2) == where_maps_to(max_i,min_j,2) .and. &
3518           where_maps_to(min_i,min_j,1) /= OUTSIDE_DOMAIN) then
3519          x_dest = where_maps_to(min_i,min_j,1)
3520          y_dest = where_maps_to(min_i,min_j,2)
3521
3522          ! If this grid point was already given a value from higher-priority source data,
3523          !   there is nothing to do.
3524 !         if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
3525
3526             ! If this grid point has never been given a value by this level of source data,
3527             !   initialize the point
3528             if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
3529                do k=min_k,max_k
3530                   dst_array(x_dest,y_dest,k) = 0.
3531                end do
3532             end if
3533
3534             ! Sum all the points whose nearest neighbor is this grid point
3535             if (present(mask_array) .and. present(mask_relational)) then
3536                do i=min_i,max_i
3537                   do j=min_j,max_j
3538                      do k=min_k,max_k
3539                         ! Ignore masked/missing values in the source data
3540                         if (tile_array(i,j,k) /= msgval) then
3541                            if (mask_relational == ' ' .and. mask_array(i,j) /= maskval) then
3542                               dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3543                               n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3544                               call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3545                            else if (mask_relational == '<' .and. mask_array(i,j) >= maskval) then
3546                               dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3547                               n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3548                               call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3549                            else if (mask_relational == '>' .and. mask_array(i,j) <= maskval) then
3550                               dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3551                               n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3552                               call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3553                            end if
3554                         end if
3555                      end do
3556                   end do
3557                end do
3558             else if (present(mask_array)) then
3559                do i=min_i,max_i
3560                   do j=min_j,max_j
3561                      do k=min_k,max_k
3562                         ! Ignore masked/missing values in the source data
3563                         if ((tile_array(i,j,k) /= msgval) .and. &
3564                             (mask_array(i,j) /= maskval)) then
3565                            dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3566                            n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3567                            call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3568                         end if
3569                      end do
3570                   end do
3571                end do
3572             else
3573                do i=min_i,max_i
3574                   do j=min_j,max_j
3575                      do k=min_k,max_k
3576                         ! Ignore masked/missing values in the source data
3577                         if ((tile_array(i,j,k) /= msgval)) then
3578                            dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3579                            n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3580                            call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3581                         end if
3582                      end do
3583                   end do
3584                end do
3585             end if
3586
3587 !         end if
3588
3589       ! Rectangle is a square of four points, and we can simply deal with each of the points
3590       else if (((max_i - min_i + 1) <= 2) .and. ((max_j - min_j + 1) <= 2)) then
3591          do i=min_i,max_i
3592             do j=min_j,max_j
3593                x_dest = where_maps_to(i,j,1)
3594                y_dest = where_maps_to(i,j,2)
3595
3596                if (x_dest /= OUTSIDE_DOMAIN) then
3597
3598 !                  if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
3599                      if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
3600                         do k=min_k,max_k
3601                            dst_array(x_dest,y_dest,k) = 0.
3602                         end do
3603                      end if
3604
3605                      if (present(mask_array) .and. present(mask_relational)) then
3606                         do k=min_k,max_k
3607                            ! Ignore masked/missing values
3608                            if (tile_array(i,j,k) /= msgval) then
3609                               if (mask_relational == ' ' .and. mask_array(i,j) /= maskval) then
3610                                  dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3611                                  n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3612                                  call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3613                               else if (mask_relational == '<' .and. mask_array(i,j) >= maskval) then
3614                                  dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3615                                  n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3616                                  call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3617                               else if (mask_relational == '>' .and. mask_array(i,j) <= maskval) then
3618                                  dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3619                                  n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3620                                  call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3621                               end if
3622                            end if
3623                         end do
3624                      else if (present(mask_array)) then
3625                         do k=min_k,max_k
3626                            ! Ignore masked/missing values
3627                            if ((tile_array(i,j,k) /= msgval) .and. &
3628                                 (mask_array(i,j) /= maskval)) then
3629                               dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3630                               n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3631                               call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3632                            end if
3633                         end do
3634                      else
3635                         do k=min_k,max_k
3636                            ! Ignore masked/missing values
3637                            if ((tile_array(i,j,k) /= msgval)) then
3638                               dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
3639                               n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
3640                               call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
3641                            end if
3642                         end do
3643                      end if
3644 !                  end if
3645
3646                end if
3647             end do
3648          end do
3649
3650       ! Not all corners map to the same grid point, and the rectangle contains more than
3651       !   four points
3652       else
3653          center_i = (max_i + min_i)/2
3654          center_j = (max_j + min_j)/2
3655
3656          ! Recursively process lower-left rectangle
3657          call process_continuous_block(tile_array, where_maps_to, &
3658                     src_min_x, src_min_y, src_min_z, &
3659                     src_max_x, src_max_y, src_max_z, &
3660                     min_i, min_j, min_k, &
3661                     center_i, center_j, max_k, &
3662                     dst_array, n, &
3663                     start_x, end_x, start_y, end_y, start_z, end_z, &
3664                     istagger, &
3665                     new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
3666
3667          if (center_i < max_i) then
3668             ! Recursively process lower-right rectangle
3669             call process_continuous_block(tile_array, where_maps_to, &
3670                        src_min_x, src_min_y, src_min_z, &
3671                        src_max_x, src_max_y, src_max_z, &
3672                        center_i+1, min_j, min_k, max_i, &
3673                        center_j, max_k, &
3674                        dst_array, n, &
3675                        start_x, end_x, start_y, &
3676                        end_y, start_z, end_z, &
3677                        istagger, &
3678                        new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
3679          end if
3680
3681          if (center_j < max_j) then
3682             ! Recursively process upper-left rectangle
3683             call process_continuous_block(tile_array, where_maps_to, &
3684                        src_min_x, src_min_y, src_min_z, &
3685                        src_max_x, src_max_y, src_max_z, &
3686                        min_i, center_j+1, min_k, center_i, &
3687                        max_j, max_k, &
3688                        dst_array, n, &
3689                        start_x, end_x, start_y, &
3690                        end_y, start_z, end_z, &
3691                        istagger, &
3692                        new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
3693          end if
3694
3695          if (center_i < max_i .and. center_j < max_j) then
3696             ! Recursively process upper-right rectangle
3697             call process_continuous_block(tile_array, where_maps_to, &
3698                        src_min_x, src_min_y, src_min_z, &
3699                        src_max_x, src_max_y, src_max_z, &
3700                        center_i+1, center_j+1, min_k, max_i, &
3701                        max_j, max_k, &
3702                        dst_array, n, &
3703                        start_x, end_x, start_y, &
3704                        end_y, start_z, end_z, &
3705                        istagger, &
3706                        new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
3707          end if
3708       end if
3709
3710    end subroutine process_continuous_block
3711
3712 end module process_domain_module