| blob | 8097d08f0ee218b42b621d767ebebf600839258a |
1 module process_domain_module
3 contains
5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
6 ! Name: process_domain
7 !
8 ! Purpose:
9 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
10 subroutine process_domain(n, extra_row, extra_col)
12 use date_pack
13 use gridinfo_module
14 use interp_option_module
15 use misc_definitions_module
16 use module_debug
17 use storage_module
19 implicit none
21 ! Arguments
22 integer, intent(in) :: n
23 logical, intent(in) :: extra_row, extra_col
25 ! Local variables
26 integer :: i, t, dyn_opt, &
27 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
28 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
29 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
30 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
31 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
32 idiff, n_times, &
33 west_east_dim, south_north_dim, bottom_top_dim, map_proj, &
34 is_water, is_lake, is_ice, is_urban, i_soilwater, &
35 grid_id, parent_id, i_parent_start, j_parent_start, &
36 i_parent_end, j_parent_end, parent_grid_ratio, sub_x, sub_y, num_land_cat, process_bdy_width
37 real :: cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
38 dom_dx, dom_dy, pole_lat, pole_lon
39 real, dimension(16) :: corner_lats, corner_lons
40 real, pointer, dimension(:,:) :: landmask
41 real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
42 logical, allocatable, dimension(:) :: got_this_field, got_const_field
43 character (len=19) :: valid_date, temp_date
44 character (len=128) :: title, mminlu
45 character (len=128), allocatable, dimension(:) :: output_flags, td_output_flags
46 character (len=128), dimension(:), pointer :: geogrid_flags
48 ! CWH Initialize local pointer variables
49 nullify(landmask)
50 nullify(xlat)
51 nullify(xlon)
52 nullify(xlat_u)
53 nullify(xlon_u)
54 nullify(xlat_v)
55 nullify(xlon_v)
56 nullify(geogrid_flags)
58 ! Compute number of times that we will process
59 call geth_idts(end_date(n), start_date(n), idiff)
60 call mprintf((idiff < 0),ERROR,'Ending date is earlier than starting date in namelist for domain %i.', i1=n)
62 n_times = idiff / interval_seconds
64 ! Check that the interval evenly divides the range of times to process
65 call mprintf((mod(idiff, interval_seconds) /= 0),WARN, &
66 'In namelist, interval_seconds does not evenly divide '// &
67 '(end_date - start_date) for domain %i. Only %i time periods '// &
68 'will be processed.', i1=n, i2=n_times)
70 ! Initialize the storage module
71 call mprintf(.true.,LOGFILE,'Initializing storage module')
72 call storage_init()
74 !
75 ! Do time-independent processing
76 !
77 call get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, map_proj, &
78 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
79 we_patch_s, we_patch_e, &
80 we_patch_stag_s, we_patch_stag_e, &
81 sn_patch_s, sn_patch_e, &
82 sn_patch_stag_s, sn_patch_stag_e, &
83 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
84 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
85 mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
86 grid_id, parent_id, &
87 i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
88 parent_grid_ratio, sub_x, sub_y, &
89 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
90 pole_lat, pole_lon, dom_dx, dom_dy, landmask, xlat, xlon, xlat_u, xlon_u, &
91 xlat_v, xlon_v, corner_lats, corner_lons, title, geogrid_flags)
94 allocate(output_flags(num_entries))
95 allocate(got_const_field(num_entries))
97 do i=1,num_entries
98 output_flags(i) = ' '
99 got_const_field(i) = .false.
100 end do
102 ! This call is to process the constant met fields (SST or SEAICE, for example)
103 ! That we process constant fields is indicated by the first argument
104 call process_single_met_time(.true., temp_date, n, extra_row, extra_col, xlat, xlon, &
105 xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
106 title, dyn_opt, &
107 west_east_dim, south_north_dim, &
108 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
109 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
110 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
111 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
112 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
113 got_const_field, &
114 map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
115 grid_id, parent_id, i_parent_start, &
116 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
117 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
118 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
119 corner_lats, corner_lons, output_flags, geogrid_flags, 0)
121 !
122 ! Begin time-dependent processing
123 !
125 allocate(td_output_flags(num_entries))
126 allocate(got_this_field (num_entries))
128 ! Loop over all times to be processed for this domain
129 do t=0,n_times
131 call geth_newdate(valid_date, trim(start_date(n)), t*interval_seconds)
132 temp_date = ' '
134 if (mod(interval_seconds,3600) == 0) then
135 write(temp_date,'(a13)') valid_date(1:10)//'_'//valid_date(12:13)
136 else if (mod(interval_seconds,60) == 0) then
137 write(temp_date,'(a16)') valid_date(1:10)//'_'//valid_date(12:16)
138 else
139 write(temp_date,'(a19)') valid_date(1:10)//'_'//valid_date(12:19)
140 end if
142 call mprintf(.true.,STDOUT, ' Processing %s', s1=trim(temp_date))
143 call mprintf(.true.,LOGFILE, 'Preparing to process output time %s', s1=temp_date)
145 do i=1,num_entries
146 td_output_flags(i) = output_flags(i)
147 got_this_field(i) = got_const_field(i)
148 end do
150 if (t > 0) then
151 process_bdy_width = process_only_bdy
152 else
153 process_bdy_width = 0
154 end if
156 call process_single_met_time(.false., temp_date, n, extra_row, extra_col, xlat, xlon, &
157 xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
158 title, dyn_opt, &
159 west_east_dim, south_north_dim, &
160 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
161 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
162 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
163 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
164 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
165 got_this_field, &
166 map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
167 grid_id, parent_id, i_parent_start, &
168 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
169 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
170 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
171 corner_lats, corner_lons, td_output_flags, geogrid_flags, process_bdy_width)
173 end do ! Loop over n_times
176 deallocate(td_output_flags)
177 deallocate(got_this_field)
179 deallocate(output_flags)
180 deallocate(got_const_field)
182 if (associated(geogrid_flags)) deallocate(geogrid_flags)
184 call storage_delete_all()
186 end subroutine process_domain
189 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
190 ! Name: get_static_fields
191 !
192 ! Purpose:
193 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
194 subroutine get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, &
195 map_proj, &
196 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
197 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
198 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
199 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
200 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
201 mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
202 grid_id, parent_id, &
203 i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
204 parent_grid_ratio, sub_x, sub_y, &
205 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
206 pole_lat, pole_lon, dom_dx, dom_dy, landmask, xlat, xlon, xlat_u, xlon_u, &
207 xlat_v, xlon_v, corner_lats, corner_lons, title, geogrid_flags)
209 use gridinfo_module
210 use input_module
211 use llxy_module
212 use parallel_module
213 use storage_module
214 use module_debug
215 use list_module
217 implicit none
219 ! Arguments
220 integer, intent(in) :: n
221 integer, intent(inout) :: dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, &
222 map_proj, &
223 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
224 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
225 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
226 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
227 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
228 is_water, is_lake, is_ice, is_urban, i_soilwater, 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, num_land_cat
231 real, pointer, dimension(:,:) :: landmask
232 real, intent(inout) :: cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
233 dom_dx, dom_dy, pole_lat, pole_lon
234 real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
235 real, dimension(16), intent(out) :: corner_lats, corner_lons
236 character (len=128), intent(inout) :: title, mminlu
237 character (len=128), dimension(:), pointer :: geogrid_flags
239 ! Local variables
240 integer :: istatus, i, j, k, sp1, ep1, sp2, ep2, sp3, ep3, &
241 lh_mult, rh_mult, bh_mult, th_mult, subx, suby
242 integer :: we_mem_subgrid_s, we_mem_subgrid_e, &
243 sn_mem_subgrid_s, sn_mem_subgrid_e
244 integer :: we_patch_subgrid_s, we_patch_subgrid_e, &
245 sn_patch_subgrid_s, sn_patch_subgrid_e
246 real, pointer, dimension(:,:,:) :: real_array
247 character (len=3) :: memorder
248 character (len=128) :: grid_type, datestr, cname, stagger, cunits, cdesc, ctemp
249 character (len=128), dimension(3) :: dimnames
250 type (fg_input) :: field
251 type (list) :: static_list, flag_list
252 type (list_item), dimension(:), pointer :: static_list_array
254 call list_init(static_list)
256 ! CWH Initialize local pointer variables
257 nullify(real_array)
259 ! Initialize the input module to read static input data for this domain
260 call mprintf(.true.,LOGFILE,'Opening static input file.')
261 call input_init(n, istatus)
262 call mprintf((istatus /= 0),ERROR, 'input_init(): Error opening input for domain %i.', i1=n)
264 ! Read global attributes from the static data input file
265 call mprintf(.true.,LOGFILE,'Reading static global attributes.')
266 call read_global_attrs(title, datestr, grid_type, dyn_opt, west_east_dim, &
267 south_north_dim, bottom_top_dim, &
268 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
269 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
270 map_proj, mminlu, num_land_cat, &
271 is_water, is_lake, is_ice, is_urban, i_soilwater, &
272 grid_id, parent_id, i_parent_start, &
273 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
274 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, &
275 truelat2, pole_lat, pole_lon, parent_grid_ratio, &
276 corner_lats, corner_lons, sub_x, sub_y)
278 we_dom_s = 1
279 sn_dom_s = 1
280 if (grid_type(1:1) == 'C') then
281 we_dom_e = west_east_dim - 1
282 sn_dom_e = south_north_dim - 1
283 else if (grid_type(1:1) == 'E') then
284 we_dom_e = west_east_dim
285 sn_dom_e = south_north_dim
286 end if
288 ! Given the full dimensions of this domain, find out the range of indices
289 ! that will be worked on by this processor. This information is given by
290 ! my_minx, my_miny, my_maxx, my_maxy
291 call parallel_get_tile_dims(west_east_dim, south_north_dim)
293 ! Must figure out patch dimensions from info in parallel module
294 if (nprocs > 1 .and. .not. do_tiled_input) then
296 we_patch_s = my_minx
297 we_patch_stag_s = my_minx
298 we_patch_e = my_maxx - 1
299 sn_patch_s = my_miny
300 sn_patch_stag_s = my_miny
301 sn_patch_e = my_maxy - 1
303 if (gridtype == 'C') then
304 if (my_x /= nproc_x - 1) then
305 we_patch_e = we_patch_e + 1
306 we_patch_stag_e = we_patch_e
307 else
308 we_patch_stag_e = we_patch_e + 1
309 end if
310 if (my_y /= nproc_y - 1) then
311 sn_patch_e = sn_patch_e + 1
312 sn_patch_stag_e = sn_patch_e
313 else
314 sn_patch_stag_e = sn_patch_e + 1
315 end if
316 else if (gridtype == 'E') then
317 we_patch_e = we_patch_e + 1
318 sn_patch_e = sn_patch_e + 1
319 we_patch_stag_e = we_patch_e
320 sn_patch_stag_e = sn_patch_e
321 end if
323 end if
325 ! Compute multipliers for halo width; these must be 0/1
326 if (my_x /= 0) then
327 lh_mult = 1
328 else
329 lh_mult = 0
330 end if
331 if (my_x /= (nproc_x-1)) then
332 rh_mult = 1
333 else
334 rh_mult = 0
335 end if
336 if (my_y /= 0) then
337 bh_mult = 1
338 else
339 bh_mult = 0
340 end if
341 if (my_y /= (nproc_y-1)) then
342 th_mult = 1
343 else
344 th_mult = 0
345 end if
347 we_mem_s = we_patch_s - HALO_WIDTH*lh_mult
348 we_mem_e = we_patch_e + HALO_WIDTH*rh_mult
349 sn_mem_s = sn_patch_s - HALO_WIDTH*bh_mult
350 sn_mem_e = sn_patch_e + HALO_WIDTH*th_mult
351 we_mem_stag_s = we_patch_stag_s - HALO_WIDTH*lh_mult
352 we_mem_stag_e = we_patch_stag_e + HALO_WIDTH*rh_mult
353 sn_mem_stag_s = sn_patch_stag_s - HALO_WIDTH*bh_mult
354 sn_mem_stag_e = sn_patch_stag_e + HALO_WIDTH*th_mult
356 ! Initialize a proj_info type for the destination grid projection. This will
357 ! primarily be used for rotating Earth-relative winds to grid-relative winds
358 call set_domain_projection(map_proj, stand_lon, truelat1, truelat2, &
359 dom_dx, dom_dy, dom_dx, dom_dy, west_east_dim, &
360 south_north_dim, real(west_east_dim)/2., &
361 real(south_north_dim)/2.,cen_lat, cen_lon, &
362 cen_lat, cen_lon)
364 ! Read static fields using the input module; we know that there are no more
365 ! fields to be read when read_next_field() returns a non-zero status.
366 istatus = 0
367 do while (istatus == 0)
368 call read_next_field(sp1, ep1, sp2, ep2, sp3, ep3, cname, cunits, cdesc, &
369 memorder, stagger, dimnames, subx, suby, &
370 real_array, istatus)
371 if (istatus == 0) then
373 call mprintf(.true.,LOGFILE, 'Read in static field %s.',s1=cname)
374 call list_insert(static_list, ckey=cname, cvalue=cname)
376 ! We will also keep copies in core of the lat/lon arrays, for use in
377 ! interpolation of the met fields to the model grid.
378 ! For now, we assume that the lat/lon arrays will have known field names
379 if (index(cname, 'XLAT_M') /= 0 .and. &
380 len_trim(cname) == len_trim('XLAT_M')) then
381 allocate(xlat(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
382 xlat(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
383 call exchange_halo_r(xlat, &
384 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
385 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
387 else if (index(cname, 'XLONG_M') /= 0 .and. &
388 len_trim(cname) == len_trim('XLONG_M')) then
389 allocate(xlon(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
390 xlon(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
391 call exchange_halo_r(xlon, &
392 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
393 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
395 else if (index(cname, 'XLAT_U') /= 0 .and. &
396 len_trim(cname) == len_trim('XLAT_U')) then
397 allocate(xlat_u(we_mem_stag_s:we_mem_stag_e,sn_mem_s:sn_mem_e))
398 xlat_u(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
399 call exchange_halo_r(xlat_u, &
400 we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
401 we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
403 else if (index(cname, 'XLONG_U') /= 0 .and. &
404 len_trim(cname) == len_trim('XLONG_U')) then
405 allocate(xlon_u(we_mem_stag_s:we_mem_stag_e,sn_mem_s:sn_mem_e))
406 xlon_u(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
407 call exchange_halo_r(xlon_u, &
408 we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
409 we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
411 else if (index(cname, 'XLAT_V') /= 0 .and. &
412 len_trim(cname) == len_trim('XLAT_V')) then
413 allocate(xlat_v(we_mem_s:we_mem_e,sn_mem_stag_s:sn_mem_stag_e))
414 xlat_v(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(:,:,1)
415 call exchange_halo_r(xlat_v, &
416 we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
417 we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
419 else if (index(cname, 'XLONG_V') /= 0 .and. &
420 len_trim(cname) == len_trim('XLONG_V')) then
421 allocate(xlon_v(we_mem_s:we_mem_e,sn_mem_stag_s:sn_mem_stag_e))
422 xlon_v(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(:,:,1)
423 call exchange_halo_r(xlon_v, &
424 we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
425 we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
427 else if (index(cname, 'LANDMASK') /= 0 .and. &
428 len_trim(cname) == len_trim('LANDMASK')) then
429 allocate(landmask(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
430 landmask(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
431 call exchange_halo_r(landmask, &
432 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
433 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
435 end if
437 if (subx > 1) then
438 we_mem_subgrid_s = (we_mem_s + HALO_WIDTH*lh_mult - 1)*subx - HALO_WIDTH*lh_mult + 1
439 we_mem_subgrid_e = (we_mem_e + (1-rh_mult) - HALO_WIDTH*rh_mult )*subx + HALO_WIDTH*rh_mult
440 we_patch_subgrid_s = (we_patch_s - 1)*subx + 1
441 we_patch_subgrid_e = (we_patch_e + (1-rh_mult) )*subx
442 end if
443 if (suby > 1) then
444 sn_mem_subgrid_s = (sn_mem_s + HALO_WIDTH*bh_mult - 1)*suby - HALO_WIDTH*bh_mult + 1
445 sn_mem_subgrid_e = (sn_mem_e + (1-th_mult) - HALO_WIDTH*th_mult )*suby + HALO_WIDTH*th_mult
446 sn_patch_subgrid_s = (sn_patch_s - 1)*suby + 1
447 sn_patch_subgrid_e = (sn_patch_e + (1-th_mult) )*suby
448 end if
450 ! Having read in a field, we write each level individually to the
451 ! storage module; levels will be reassembled later on when they
452 ! are written.
453 do k=sp3,ep3
454 field%header%sr_x=subx
455 field%header%sr_y=suby
456 field%header%version = 1
457 field%header%date = start_date(n)
458 field%header%time_dependent = .false.
459 field%header%mask_field = .false.
460 field%header%constant_field = .false.
461 field%header%forecast_hour = 0.0
462 field%header%fg_source = 'geogrid_model'
463 field%header%field = cname
464 field%header%units = cunits
465 field%header%description = cdesc
466 field%header%vertical_coord = dimnames(3)
467 field%header%vertical_level = k
468 field%header%array_order = memorder
469 field%header%is_wind_grid_rel = .true.
470 field%header%array_has_missing_values = .false.
471 if (gridtype == 'C') then
472 if (subx > 1 .or. suby > 1) then
473 field%map%stagger = M
474 field%header%dim1(1) = we_mem_subgrid_s
475 field%header%dim1(2) = we_mem_subgrid_e
476 field%header%dim2(1) = sn_mem_subgrid_s
477 field%header%dim2(2) = sn_mem_subgrid_e
478 else if (trim(stagger) == 'M') then
479 field%map%stagger = M
480 field%header%dim1(1) = we_mem_s
481 field%header%dim1(2) = we_mem_e
482 field%header%dim2(1) = sn_mem_s
483 field%header%dim2(2) = sn_mem_e
484 else if (trim(stagger) == 'U') then
485 field%map%stagger = U
486 field%header%dim1(1) = we_mem_stag_s
487 field%header%dim1(2) = we_mem_stag_e
488 field%header%dim2(1) = sn_mem_s
489 field%header%dim2(2) = sn_mem_e
490 else if (trim(stagger) == 'V') then
491 field%map%stagger = V
492 field%header%dim1(1) = we_mem_s
493 field%header%dim1(2) = we_mem_e
494 field%header%dim2(1) = sn_mem_stag_s
495 field%header%dim2(2) = sn_mem_stag_e
496 end if
497 else if (gridtype == 'E') then
498 if (trim(stagger) == 'M') then
499 field%map%stagger = HH
500 else if (trim(stagger) == 'V') then
501 field%map%stagger = VV
502 end if
503 field%header%dim1(1) = we_mem_s
504 field%header%dim1(2) = we_mem_e
505 field%header%dim2(1) = sn_mem_s
506 field%header%dim2(2) = sn_mem_e
507 end if
509 allocate(field%valid_mask)
511 if (subx > 1 .or. suby > 1) then
512 allocate(field%r_arr(we_mem_subgrid_s:we_mem_subgrid_e,&
513 sn_mem_subgrid_s:sn_mem_subgrid_e))
514 field%r_arr(we_patch_subgrid_s:we_patch_subgrid_e,sn_patch_subgrid_s:sn_patch_subgrid_e) = &
515 real_array(sp1:ep1,sp2:ep2,k)
516 call exchange_halo_r(field%r_arr, &
517 we_mem_subgrid_s, we_mem_subgrid_e, sn_mem_subgrid_s, sn_mem_subgrid_e, 1, 1, &
518 we_patch_subgrid_s, we_patch_subgrid_e, sn_patch_subgrid_s, sn_patch_subgrid_e, 1, 1)
519 call bitarray_create(field%valid_mask, &
520 (we_mem_subgrid_e-we_mem_subgrid_s)+1, &
521 (sn_mem_subgrid_e-sn_mem_subgrid_s)+1)
522 do j=1,(sn_mem_subgrid_e-sn_mem_subgrid_s)+1
523 do i=1,(we_mem_subgrid_e-we_mem_subgrid_s)+1
524 call bitarray_set(field%valid_mask, i, j)
525 end do
526 end do
528 else if (field%map%stagger == M .or. &
529 field%map%stagger == HH .or. &
530 field%map%stagger == VV) then
531 allocate(field%r_arr(we_mem_s:we_mem_e,&
532 sn_mem_s:sn_mem_e))
533 field%r_arr(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(sp1:ep1,sp2:ep2,k)
534 call exchange_halo_r(field%r_arr, &
535 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
536 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
537 call bitarray_create(field%valid_mask, &
538 (we_mem_e-we_mem_s)+1, &
539 (sn_mem_e-sn_mem_s)+1)
540 do j=1,(sn_mem_e-sn_mem_s)+1
541 do i=1,(we_mem_e-we_mem_s)+1
542 call bitarray_set(field%valid_mask, i, j)
543 end do
544 end do
545 else if (field%map%stagger == U) then
546 allocate(field%r_arr(we_mem_stag_s:we_mem_stag_e,&
547 sn_mem_s:sn_mem_e))
548 field%r_arr(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(sp1:ep1,sp2:ep2,k)
549 call exchange_halo_r(field%r_arr, &
550 we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
551 we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
552 call bitarray_create(field%valid_mask, &
553 (we_mem_stag_e-we_mem_stag_s)+1, &
554 (sn_mem_e-sn_mem_s)+1)
555 do j=1,(sn_mem_e-sn_mem_s)+1
556 do i=1,(we_mem_stag_e-we_mem_stag_s)+1
557 call bitarray_set(field%valid_mask, i, j)
558 end do
559 end do
560 else if (field%map%stagger == V) then
561 allocate(field%r_arr(we_mem_s:we_mem_e,&
562 sn_mem_stag_s:sn_mem_stag_e))
563 field%r_arr(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(sp1:ep1,sp2:ep2,k)
564 call exchange_halo_r(field%r_arr, &
565 we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
566 we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
567 call bitarray_create(field%valid_mask, &
568 (we_mem_e-we_mem_s)+1, &
569 (sn_mem_stag_e-sn_mem_stag_s)+1)
570 do j=1,(sn_mem_stag_e-sn_mem_stag_s)+1
571 do i=1,(we_mem_e-we_mem_s)+1
572 call bitarray_set(field%valid_mask, i, j)
573 end do
574 end do
575 end if
577 nullify(field%modified_mask)
579 call storage_put_field(field)
581 end do
583 end if
584 end do
586 static_list_array => list_get_keys(static_list)
587 call list_init(flag_list)
588 do i=1,size(static_list_array)
589 istatus = 0
590 ctemp = 'FLAG_'//trim(static_list_array(i)%ckey)
591 call ext_get_dom_ti_integer_scalar(ctemp, istatus, suppress_errors=.true.)
592 if (istatus == 1) call list_insert(flag_list, ckey=ctemp, cvalue=ctemp)
593 end do
594 deallocate(static_list_array)
595 call list_destroy(static_list)
597 ! Done reading all static fields for this domain
598 call input_close()
600 static_list_array => list_get_keys(flag_list)
601 allocate(geogrid_flags(size(static_list_array)))
602 do i=1,size(static_list_array)
603 geogrid_flags(i) = static_list_array(i)%ckey
604 end do
605 deallocate(static_list_array)
607 call list_destroy(flag_list)
609 end subroutine get_static_fields
612 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
613 ! Name: process_single_met_time
614 !
615 ! Purpose:
616 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
617 subroutine process_single_met_time(do_const_processing, &
618 temp_date, n, extra_row, extra_col, xlat, xlon, &
619 xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
620 title, dyn_opt, &
621 west_east_dim, south_north_dim, &
622 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
623 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
624 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
625 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
626 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
627 got_this_field, &
628 map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
629 grid_id, parent_id, i_parent_start, &
630 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
631 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
632 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
633 corner_lats, corner_lons, output_flags, geogrid_flags, process_bdy_width)
635 use bitarray_module
636 use gridinfo_module
637 use interp_module
638 use interp_option_module
639 use llxy_module
640 use misc_definitions_module
641 use module_debug
642 use output_module
643 use parallel_module
644 use read_met_module
645 use rotate_winds_module
646 use storage_module
648 implicit none
650 ! Arguments
651 logical, intent(in) :: do_const_processing
652 integer, intent(in) :: n, dyn_opt, west_east_dim, south_north_dim, map_proj, &
653 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
654 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
655 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
656 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
657 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
658 is_water, is_lake, is_ice, is_urban, i_soilwater, &
659 grid_id, parent_id, i_parent_start, j_parent_start, &
660 i_parent_end, j_parent_end, parent_grid_ratio, sub_x, sub_y, num_land_cat, &
661 process_bdy_width
662 ! BUG: Should we be passing these around as pointers, or just declare them as arrays?
663 real, pointer, dimension(:,:) :: landmask
664 real, intent(in) :: dom_dx, dom_dy, cen_lat, moad_cen_lat, cen_lon, stand_lon, &
665 truelat1, truelat2, pole_lat, pole_lon
666 real, dimension(16), intent(in) :: corner_lats, corner_lons
667 real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
668 logical, intent(in) :: extra_row, extra_col
669 logical, dimension(:), intent(inout) :: got_this_field
670 character (len=19), intent(in) :: temp_date
671 character (len=128), intent(in) :: mminlu
672 character (len=128), dimension(:), intent(inout) :: output_flags
673 character (len=128), dimension(:), pointer :: geogrid_flags
675 ! BUG: Move this constant to misc_definitions_module?
676 integer, parameter :: BDR_WIDTH = 3
678 ! Local variables
679 integer :: istatus, iqstatus, fg_idx, idx, idxt, i, j, bottom_top_dim, &
680 sm1, em1, sm2, em2, sm3, em3, &
681 sp1, ep1, sp2, ep2, sp3, ep3, &
682 sd1, ed1, sd2, ed2, sd3, ed3, &
683 u_idx, bdr_wdth
684 integer :: nmet_flags
685 integer :: num_metgrid_soil_levs
686 integer, pointer, dimension(:) :: soil_levels
687 real :: rx, ry
688 real :: threshold
689 logical :: do_gcell_interp
690 integer, pointer, dimension(:) :: u_levels, v_levels
691 real, pointer, dimension(:,:) :: halo_slab
692 real, pointer, dimension(:,:,:) :: real_array
693 character (len=19) :: output_date
694 character (len=128) :: cname, title
695 character (len=MAX_FILENAME_LEN) :: input_name
696 character (len=128), allocatable, dimension(:) :: met_flags
697 type (fg_input) :: field, u_field, v_field
698 type (met_data) :: fg_data
700 ! CWH Initialize local pointer variables
701 nullify(soil_levels)
702 nullify(u_levels)
703 nullify(v_levels)
704 nullify(halo_slab)
705 nullify(real_array)
708 ! For this time, we need to process all first-guess filename roots. When we
709 ! hit a root containing a '*', we assume we have hit the end of the list
710 fg_idx = 1
711 if (do_const_processing) then
712 input_name = constants_name(fg_idx)
713 else
714 input_name = fg_name(fg_idx)
715 end if
716 do while (input_name /= '*')
718 call mprintf(.true.,STDOUT, ' %s', s1=input_name)
719 call mprintf(.true.,LOGFILE, 'Getting input fields from %s', s1=input_name)
721 ! Do a first pass through this fg source to get all mask fields used
722 ! during interpolation
723 call get_interp_masks(trim(input_name), do_const_processing, temp_date)
725 istatus = 0
727 ! Initialize the module for reading in the met fields
728 call read_met_init(trim(input_name), do_const_processing, temp_date, istatus)
730 if (istatus == 0) then
732 ! Process all fields and levels from the current file; read_next_met_field()
733 ! will return a non-zero status when there are no more fields to be read.
734 do while (istatus == 0)
736 call read_next_met_field(fg_data, istatus)
738 if (istatus == 0) then
740 ! Find index into fieldname, interp_method, masked, and fill_missing
741 ! of the current field
742 idxt = num_entries + 1
743 do idx=1,num_entries
744 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
745 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
747 got_this_field(idx) = .true.
749 if (index(input_name,trim(from_input(idx))) /= 0 .or. &
750 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
751 idxt = idx
752 end if
754 end if
755 end do
756 idx = idxt
757 if (idx > num_entries) idx = num_entries ! The last entry is a default
759 ! Do we need to rename this field?
760 if (output_name(idx) /= ' ') then
761 fg_data%field = output_name(idx)(1:9)
763 idxt = num_entries + 1
764 do idx=1,num_entries
765 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
766 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
768 got_this_field(idx) = .true.
770 if (index(input_name,trim(from_input(idx))) /= 0 .or. &
771 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
772 idxt = idx
773 end if
775 end if
776 end do
777 idx = idxt
778 if (idx > num_entries) idx = num_entries ! The last entry is a default
779 end if
781 ! Do a simple check to see whether this is a global dataset
782 ! Note that we do not currently support regional Gaussian grids
783 if ((fg_data%iproj == PROJ_LATLON .and. abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) &
784 .or. (fg_data%iproj == PROJ_GAUSS)) then
785 bdr_wdth = BDR_WIDTH
786 allocate(halo_slab(1-BDR_WIDTH:fg_data%nx+BDR_WIDTH,1:fg_data%ny))
788 halo_slab(1:fg_data%nx, 1:fg_data%ny) = &
789 fg_data%slab(1:fg_data%nx, 1:fg_data%ny)
791 halo_slab(1-BDR_WIDTH:0, 1:fg_data%ny) = &
792 fg_data%slab(fg_data%nx-BDR_WIDTH+1:fg_data%nx, 1:fg_data%ny)
794 halo_slab(fg_data%nx+1:fg_data%nx+BDR_WIDTH, 1:fg_data%ny) = &
795 fg_data%slab(1:BDR_WIDTH, 1:fg_data%ny)
797 deallocate(fg_data%slab)
798 else
799 bdr_wdth = 0
800 halo_slab => fg_data%slab
801 nullify(fg_data%slab)
802 end if
804 call mprintf(.true.,LOGFILE,'Processing %s at level %f.',s1=fg_data%field,f1=fg_data%xlvl)
806 call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, &
807 fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, &
808 fg_data%deltalon, fg_data%starti, fg_data%startj, &
809 fg_data%startlat, fg_data%startlon, earth_radius=fg_data%earth_radius*1000.)
811 ! Initialize fg_input structure to store the field
812 field%header%version = 1
813 field%header%date = fg_data%hdate//' '
814 if (do_const_processing) then
815 field%header%time_dependent = .false.
816 field%header%constant_field = .true.
817 else
818 field%header%time_dependent = .true.
819 field%header%constant_field = .false.
820 end if
821 field%header%forecast_hour = fg_data%xfcst
822 field%header%fg_source = 'FG'
823 field%header%field = ' '
824 field%header%field(1:9) = fg_data%field
825 field%header%units = ' '
826 field%header%units(1:25) = fg_data%units
827 field%header%description = ' '
828 field%header%description(1:46) = fg_data%desc
829 call get_z_dim_name(fg_data%field,field%header%vertical_coord)
830 field%header%vertical_level = nint(fg_data%xlvl)
831 field%header%sr_x = 1
832 field%header%sr_y = 1
833 field%header%array_order = 'XY '
834 field%header%is_wind_grid_rel = fg_data%is_wind_grid_rel
835 field%header%array_has_missing_values = .false.
836 nullify(field%r_arr)
837 nullify(field%valid_mask)
838 nullify(field%modified_mask)
840 if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
841 output_flags(idx) = flag_in_output(idx)
842 end if
844 ! If we should not output this field, just list it as a mask field
845 if (output_this_field(idx)) then
846 field%header%mask_field = .false.
847 else
848 field%header%mask_field = .true.
849 end if
851 !
852 ! Before actually doing any interpolation to the model grid, we must check
853 ! whether we will be using the average_gcell interpolator that averages all
854 ! source points in each model grid cell
855 !
856 do_gcell_interp = .false.
857 if (index(interp_method(idx),'average_gcell') /= 0) then
859 call get_gcell_threshold(interp_method(idx), threshold, istatus)
860 if (istatus == 0) then
861 if (fg_data%dx == 0. .and. fg_data%dy == 0. .and. &
862 fg_data%deltalat /= 0. .and. fg_data%deltalon /= 0.) then
863 fg_data%dx = abs(fg_data%deltalon)
864 fg_data%dy = abs(fg_data%deltalat)
865 else
866 ! BUG: Need to more correctly handle dx/dy in meters.
867 fg_data%dx = fg_data%dx / 111000. ! Convert meters to approximate degrees
868 fg_data%dy = fg_data%dy / 111000.
869 end if
870 if (gridtype == 'C') then
871 if (threshold*max(fg_data%dx,fg_data%dy)*111. <= max(dom_dx,dom_dy)/1000.) &
872 do_gcell_interp = .true.
873 else if (gridtype == 'E') then
874 if (threshold*max(fg_data%dx,fg_data%dy) <= max(dom_dx,dom_dy)) &
875 do_gcell_interp = .true.
876 end if
877 end if
878 end if
880 ! Interpolate to U staggering
881 if (output_stagger(idx) == U) then
883 call storage_query_field(field, iqstatus)
884 if (iqstatus == 0) then
885 call storage_get_field(field, iqstatus)
886 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
887 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
888 if (associated(field%modified_mask)) then
889 call bitarray_destroy(field%modified_mask)
890 nullify(field%modified_mask)
891 end if
892 else
893 allocate(field%valid_mask)
894 call bitarray_create(field%valid_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
895 end if
897 ! Save a copy of the fg_input structure for the U field so that we can find it later
898 if (is_u_field(idx)) call dup(field, u_field)
900 allocate(field%modified_mask)
901 call bitarray_create(field%modified_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
903 if (do_const_processing .or. field%header%time_dependent) then
904 call interp_met_field(input_name, fg_data%field, U, M, &
905 field, xlat_u, xlon_u, we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, &
906 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
907 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
908 field%modified_mask, process_bdy_width)
909 else
910 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
911 end if
913 ! Interpolate to V staggering
914 else if (output_stagger(idx) == V) then
916 call storage_query_field(field, iqstatus)
917 if (iqstatus == 0) then
918 call storage_get_field(field, iqstatus)
919 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
920 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
921 if (associated(field%modified_mask)) then
922 call bitarray_destroy(field%modified_mask)
923 nullify(field%modified_mask)
924 end if
925 else
926 allocate(field%valid_mask)
927 call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
928 end if
930 ! Save a copy of the fg_input structure for the V field so that we can find it later
931 if (is_v_field(idx)) call dup(field, v_field)
933 allocate(field%modified_mask)
934 call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
936 if (do_const_processing .or. field%header%time_dependent) then
937 call interp_met_field(input_name, fg_data%field, V, M, &
938 field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
939 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
940 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
941 field%modified_mask, process_bdy_width)
942 else
943 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
944 end if
946 ! Interpolate to VV staggering
947 else if (output_stagger(idx) == VV) then
949 call storage_query_field(field, iqstatus)
950 if (iqstatus == 0) then
951 call storage_get_field(field, iqstatus)
952 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
953 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
954 if (associated(field%modified_mask)) then
955 call bitarray_destroy(field%modified_mask)
956 nullify(field%modified_mask)
957 end if
958 else
959 allocate(field%valid_mask)
960 call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
961 end if
963 ! Save a copy of the fg_input structure for the U field so that we can find it later
964 if (is_u_field(idx)) call dup(field, u_field)
966 ! Save a copy of the fg_input structure for the V field so that we can find it later
967 if (is_v_field(idx)) call dup(field, v_field)
969 allocate(field%modified_mask)
970 call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
972 if (do_const_processing .or. field%header%time_dependent) then
973 call interp_met_field(input_name, fg_data%field, VV, M, &
974 field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
975 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
976 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
977 field%modified_mask, process_bdy_width)
978 else
979 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
980 end if
982 ! All other fields interpolated to M staggering for C grid, H staggering for E grid
983 else
985 call storage_query_field(field, iqstatus)
986 if (iqstatus == 0) then
987 call storage_get_field(field, iqstatus)
988 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
989 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
990 if (associated(field%modified_mask)) then
991 call bitarray_destroy(field%modified_mask)
992 nullify(field%modified_mask)
993 end if
994 else
995 allocate(field%valid_mask)
996 call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
997 end if
999 allocate(field%modified_mask)
1000 call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
1002 if (do_const_processing .or. field%header%time_dependent) then
1003 if (gridtype == 'C') then
1004 call interp_met_field(input_name, fg_data%field, M, M, &
1005 field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1006 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1007 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
1008 field%modified_mask, process_bdy_width, landmask)
1010 else if (gridtype == 'E') then
1011 call interp_met_field(input_name, fg_data%field, HH, M, &
1012 field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1013 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1014 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
1015 field%modified_mask, process_bdy_width, landmask)
1016 end if
1017 else
1018 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
1019 end if
1021 end if
1023 call bitarray_merge(field%valid_mask, field%modified_mask)
1025 deallocate(halo_slab)
1027 ! Store the interpolated field
1028 call storage_put_field(field)
1030 call pop_source_projection()
1032 end if
1033 end do
1035 call read_met_close()
1037 call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, &
1038 fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, &
1039 fg_data%deltalon, fg_data%starti, fg_data%startj, &
1040 fg_data%startlat, fg_data%startlon, earth_radius=fg_data%earth_radius*1000.)
1042 !
1043 ! If necessary, rotate winds to earth-relative for this fg source
1044 !
1046 call storage_get_levels(u_field, u_levels)
1047 call storage_get_levels(v_field, v_levels)
1049 if (associated(u_levels) .and. associated(v_levels)) then
1050 u_idx = 1
1051 do u_idx = 1, size(u_levels)
1052 u_field%header%vertical_level = u_levels(u_idx)
1053 call storage_get_field(u_field, istatus)
1054 v_field%header%vertical_level = v_levels(u_idx)
1055 call storage_get_field(v_field, istatus)
1057 if (associated(u_field%modified_mask) .and. &
1058 associated(v_field%modified_mask)) then
1060 if (u_field%header%is_wind_grid_rel) then
1061 if (gridtype == 'C') then
1062 call map_to_met(u_field%r_arr, u_field%modified_mask, &
1063 v_field%r_arr, v_field%modified_mask, &
1064 we_mem_stag_s, sn_mem_s, &
1065 we_mem_stag_e, sn_mem_e, &
1066 we_mem_s, sn_mem_stag_s, &
1067 we_mem_e, sn_mem_stag_e, &
1068 xlon_u, xlon_v, xlat_u, xlat_v)
1069 else if (gridtype == 'E') then
1070 call map_to_met_nmm(u_field%r_arr, u_field%modified_mask, &
1071 v_field%r_arr, v_field%modified_mask, &
1072 we_mem_s, sn_mem_s, &
1073 we_mem_e, sn_mem_e, &
1074 xlat_v, xlon_v)
1075 end if
1076 end if
1078 call bitarray_destroy(u_field%modified_mask)
1079 call bitarray_destroy(v_field%modified_mask)
1080 nullify(u_field%modified_mask)
1081 nullify(v_field%modified_mask)
1082 call storage_put_field(u_field)
1083 call storage_put_field(v_field)
1084 end if
1086 end do
1088 deallocate(u_levels)
1089 deallocate(v_levels)
1091 end if
1093 call pop_source_projection()
1095 else
1096 if (do_const_processing) then
1097 call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=input_name)
1098 else
1099 call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=trim(input_name)//':'//trim(temp_date))
1100 end if
1101 end if
1103 fg_idx = fg_idx + 1
1104 if (do_const_processing) then
1105 input_name = constants_name(fg_idx)
1106 else
1107 input_name = fg_name(fg_idx)
1108 end if
1109 end do ! while (input_name /= '*')
1111 !
1112 ! Rotate winds from earth-relative to grid-relative
1113 !
1115 call storage_get_levels(u_field, u_levels)
1116 call storage_get_levels(v_field, v_levels)
1118 if (associated(u_levels) .and. associated(v_levels)) then
1119 u_idx = 1
1120 do u_idx = 1, size(u_levels)
1121 u_field%header%vertical_level = u_levels(u_idx)
1122 call storage_get_field(u_field, istatus)
1123 v_field%header%vertical_level = v_levels(u_idx)
1124 call storage_get_field(v_field, istatus)
1126 if (gridtype == 'C') then
1127 call met_to_map(u_field%r_arr, u_field%valid_mask, &
1128 v_field%r_arr, v_field%valid_mask, &
1129 we_mem_stag_s, sn_mem_s, &
1130 we_mem_stag_e, sn_mem_e, &
1131 we_mem_s, sn_mem_stag_s, &
1132 we_mem_e, sn_mem_stag_e, &
1133 xlon_u, xlon_v, xlat_u, xlat_v)
1134 else if (gridtype == 'E') then
1135 call met_to_map_nmm(u_field%r_arr, u_field%valid_mask, &
1136 v_field%r_arr, v_field%valid_mask, &
1137 we_mem_s, sn_mem_s, &
1138 we_mem_e, sn_mem_e, &
1139 xlat_v, xlon_v)
1140 end if
1142 end do
1144 deallocate(u_levels)
1145 deallocate(v_levels)
1147 end if
1149 if (do_const_processing) return
1151 !
1152 ! Now that we have all degribbed fields, we build a 3-d pressure field, and fill in any
1153 ! missing levels in the other 3-d fields
1154 !
1155 call mprintf(.true.,LOGFILE,'Filling missing levels.')
1156 call fill_missing_levels(output_flags)
1158 call mprintf(.true.,LOGFILE,'Creating derived fields.')
1159 call create_derived_fields(gridtype, fg_data%hdate, fg_data%xfcst, &
1160 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1161 we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e, &
1162 got_this_field, output_flags)
1164 !
1165 ! Check that every mandatory field was found in input data
1166 !
1167 do i=1,num_entries
1168 if (is_mandatory(i) .and. .not. got_this_field(i)) then
1169 call mprintf(.true.,ERROR,'The mandatory field %s was not found in any input data.',s1=fieldname(i))
1170 end if
1171 end do
1173 !
1174 ! Before we begin to write fields, if debug_level is set high enough, we
1175 ! write a table of which fields are available at which levels to the
1176 ! metgrid.log file, and then we check to see if any fields are not
1177 ! completely covered with data.
1178 !
1179 call storage_print_fields()
1180 call find_missing_values()
1182 !
1183 ! All of the processing is now done for this time period for this domain;
1184 ! now we simply output every field from the storage module.
1185 !
1187 title = 'OUTPUT FROM METGRID V3.6.1'
1189 ! Initialize the output module for this domain and time
1190 call mprintf(.true.,LOGFILE,'Initializing output module.')
1191 output_date = temp_date
1192 if ( .not. nocolons ) then
1193 if (len_trim(temp_date) == 13) then
1194 output_date(14:19) = ':00:00'
1195 else if (len_trim(temp_date) == 16) then
1196 output_date(17:19) = ':00'
1197 end if
1198 else
1199 if (len_trim(temp_date) == 13) then
1200 output_date(14:19) = '_00_00'
1201 else if (len_trim(temp_date) == 16) then
1202 output_date(17:19) = '_00'
1203 end if
1204 endif
1206 call output_init(n, title, output_date, gridtype, dyn_opt, &
1207 corner_lats, corner_lons, &
1208 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1209 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, &
1210 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1211 extra_col, extra_row)
1213 call get_bottom_top_dim(bottom_top_dim)
1215 ! Add in a flag to tell real that we have seven new msf fields
1216 nmet_flags = num_entries + 1
1217 if (associated(geogrid_flags)) nmet_flags = nmet_flags + size(geogrid_flags)
1218 allocate(met_flags(nmet_flags))
1219 do i=1,num_entries
1220 met_flags(i) = output_flags(i)
1221 end do
1222 if (gridtype == 'C') then
1223 met_flags(num_entries+1) = 'FLAG_MF_XY'
1224 else
1225 met_flags(num_entries+1) = ' '
1226 end if
1227 if (associated(geogrid_flags)) then
1228 do i=1,size(geogrid_flags)
1229 met_flags(num_entries+1+i) = geogrid_flags(i)
1230 end do
1231 end if
1233 ! Find out how many soil levels or layers we have; this assumes a field named ST
1234 field % header % field = 'ST'
1235 nullify(soil_levels)
1236 call storage_get_levels(field, soil_levels)
1238 if (.not. associated(soil_levels)) then
1239 field % header % field = 'SOILT'
1240 nullify(soil_levels)
1241 call storage_get_levels(field, soil_levels)
1242 end if
1244 if (.not. associated(soil_levels)) then
1245 field % header % field = 'STC_WPS'
1246 nullify(soil_levels)
1247 call storage_get_levels(field, soil_levels)
1248 end if
1250 if (associated(soil_levels)) then
1251 num_metgrid_soil_levs = size(soil_levels)
1252 deallocate(soil_levels)
1253 else
1254 num_metgrid_soil_levs = 0
1255 end if
1257 ! First write out global attributes
1258 call mprintf(.true.,LOGFILE,'Writing global attributes to output.')
1259 call write_global_attrs(title, output_date, gridtype, dyn_opt, west_east_dim, &
1260 south_north_dim, bottom_top_dim, &
1261 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
1262 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
1263 map_proj, mminlu, num_land_cat, &
1264 is_water, is_lake, is_ice, is_urban, i_soilwater, &
1265 grid_id, parent_id, i_parent_start, &
1266 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
1267 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
1268 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
1269 corner_lats, corner_lons, num_metgrid_soil_levs, &
1270 met_flags, nmet_flags, process_bdy_width)
1272 deallocate(met_flags)
1274 call reset_next_field()
1276 istatus = 0
1278 ! Now loop over all output fields, writing each to the output module
1279 do while (istatus == 0)
1280 call get_next_output_field(cname, real_array, &
1281 sm1, em1, sm2, em2, sm3, em3, istatus)
1282 if (istatus == 0) then
1284 call mprintf(.true.,LOGFILE,'Writing field %s to output.',s1=cname)
1285 call write_field(sm1, em1, sm2, em2, sm3, em3, &
1286 cname, output_date, real_array)
1287 deallocate(real_array)
1289 end if
1291 end do
1293 call mprintf(.true.,LOGFILE,'Closing output file.')
1294 call output_close()
1296 ! Free up memory used by met fields for this valid time
1297 call storage_delete_all_td()
1299 end subroutine process_single_met_time
1302 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1303 ! Name: get_interp_masks
1304 !
1305 ! Purpose:
1306 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1307 subroutine get_interp_masks(fg_prefix, is_constants, fg_date)
1309 use interp_option_module
1310 use read_met_module
1311 use storage_module
1313 implicit none
1315 ! Arguments
1316 logical, intent(in) :: is_constants
1317 character (len=*), intent(in) :: fg_prefix, fg_date
1319 ! BUG: Move this constant to misc_definitions_module?
1320 integer, parameter :: BDR_WIDTH = 3
1322 ! Local variables
1323 integer :: i, istatus, idx, idxt
1324 type (fg_input) :: mask_field
1325 type (met_data) :: fg_data
1327 istatus = 0
1329 call read_met_init(fg_prefix, is_constants, fg_date, istatus)
1331 do while (istatus == 0)
1333 call read_next_met_field(fg_data, istatus)
1335 if (istatus == 0) then
1337 ! Find out which METGRID.TBL entry goes with this field
1338 idxt = num_entries + 1
1339 do idx=1,num_entries
1340 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
1341 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
1343 if (index(fg_prefix,trim(from_input(idx))) /= 0 .or. &
1344 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1345 idxt = idx
1346 end if
1348 end if
1349 end do
1350 idx = idxt
1351 if (idx > num_entries) idx = num_entries ! The last entry is a default
1353 ! Do we need to rename this field?
1354 if (output_name(idx) /= ' ') then
1355 fg_data%field = output_name(idx)(1:9)
1357 idxt = num_entries + 1
1358 do idx=1,num_entries
1359 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
1360 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
1362 if (index(fg_prefix,trim(from_input(idx))) /= 0 .or. &
1363 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1364 idxt = idx
1365 end if
1367 end if
1368 end do
1369 idx = idxt
1370 if (idx > num_entries) idx = num_entries ! The last entry is a default
1371 end if
1373 do i=1,num_entries
1374 if (interp_mask(i) /= ' ' .and. (trim(interp_mask(i)) == trim(fg_data%field))) then
1376 mask_field%header%version = 1
1377 mask_field%header%date = ' '
1378 mask_field%header%date = fg_date
1379 if (is_constants) then
1380 mask_field%header%time_dependent = .false.
1381 mask_field%header%constant_field = .true.
1382 else
1383 mask_field%header%time_dependent = .true.
1384 mask_field%header%constant_field = .false.
1385 end if
1386 mask_field%header%mask_field = .true.
1387 mask_field%header%forecast_hour = 0.
1388 mask_field%header%fg_source = 'degribbed met data'
1389 mask_field%header%field = trim(fg_data%field)//'.mask'
1390 mask_field%header%units = '-'
1391 mask_field%header%description = '-'
1392 mask_field%header%vertical_coord = 'none'
1393 mask_field%header%vertical_level = 1
1394 mask_field%header%sr_x = 1
1395 mask_field%header%sr_y = 1
1396 mask_field%header%array_order = 'XY'
1397 mask_field%header%dim1(1) = 1
1398 mask_field%header%dim1(2) = fg_data%nx
1399 mask_field%header%dim2(1) = 1
1400 mask_field%header%dim2(2) = fg_data%ny
1401 mask_field%header%is_wind_grid_rel = .true.
1402 mask_field%header%array_has_missing_values = .false.
1403 mask_field%map%stagger = M
1405 ! Do a simple check to see whether this is a global lat/lon dataset
1406 ! Note that we do not currently support regional Gaussian grids
1407 if ((fg_data%iproj == PROJ_LATLON .and. abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) &
1408 .or. (fg_data%iproj == PROJ_GAUSS)) then
1409 allocate(mask_field%r_arr(1-BDR_WIDTH:fg_data%nx+BDR_WIDTH,1:fg_data%ny))
1411 mask_field%r_arr(1:fg_data%nx, 1:fg_data%ny) = &
1412 fg_data%slab(1:fg_data%nx, 1:fg_data%ny)
1414 mask_field%r_arr(1-BDR_WIDTH:0, 1:fg_data%ny) = &
1415 fg_data%slab(fg_data%nx-BDR_WIDTH+1:fg_data%nx, 1:fg_data%ny)
1417 mask_field%r_arr(fg_data%nx+1:fg_data%nx+BDR_WIDTH, 1:fg_data%ny) = &
1418 fg_data%slab(1:BDR_WIDTH, 1:fg_data%ny)
1420 else
1421 allocate(mask_field%r_arr(1:fg_data%nx,1:fg_data%ny))
1422 mask_field%r_arr = fg_data%slab
1423 end if
1425 nullify(mask_field%valid_mask)
1426 nullify(mask_field%modified_mask)
1428 call storage_put_field(mask_field)
1430 exit
1432 end if
1433 end do
1435 if (associated(fg_data%slab)) deallocate(fg_data%slab)
1437 end if
1439 end do
1441 call read_met_close()
1443 end subroutine get_interp_masks
1446 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1447 ! Name: interp_met_field
1448 !
1449 ! Purpose:
1450 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1451 subroutine interp_met_field(input_name, short_fieldnm, ifieldstagger, istagger, &
1452 field, xlat, xlon, sm1, em1, sm2, em2, &
1453 sd1, ed1, sd2, ed2, &
1454 slab, minx, maxx, miny, maxy, bdr, do_gcell_interp, &
1455 new_pts, process_bdy_width, landmask)
1457 use bitarray_module
1458 use interp_module
1459 use interp_option_module
1460 use llxy_module
1461 use misc_definitions_module
1462 use storage_module
1464 implicit none
1466 ! Arguments
1467 integer, intent(in) :: ifieldstagger, istagger, &
1468 sm1, em1, sm2, em2, &
1469 sd1, ed1, sd2, ed2, &
1470 minx, maxx, miny, maxy, bdr, &
1471 process_bdy_width
1472 real, dimension(minx:maxx,miny:maxy), intent(in) :: slab
1473 real, dimension(sm1:em1,sm2:em2), intent(in) :: xlat, xlon
1474 real, dimension(sm1:em1,sm2:em2), intent(in), optional :: landmask
1475 logical, intent(in) :: do_gcell_interp
1476 character (len=9), intent(in) :: short_fieldnm
1477 character (len=MAX_FILENAME_LEN), intent(in) :: input_name
1478 type (fg_input), intent(inout) :: field
1479 type (bitarray), intent(inout) :: new_pts
1481 ! Local variables
1482 integer :: i, j, idx, idxt, orig_selected_proj, interp_mask_status, &
1483 interp_land_mask_status, interp_water_mask_status, process_width
1484 integer, pointer, dimension(:) :: interp_array, interp_opts
1485 real :: rx, ry, temp
1486 real, pointer, dimension(:,:) :: data_count
1487 type (fg_input) :: mask_field, mask_water_field, mask_land_field
1488 !BPR BEGIN
1489 real, dimension(sm1:em1,sm2:em2) :: r_arr_cur_source
1490 !BPR END
1492 ! CWH Initialize local pointer variables
1493 nullify(interp_array)
1494 nullify(interp_opts)
1495 nullify(data_count)
1497 ! Find index into fieldname, interp_method, masked, and fill_missing
1498 ! of the current field
1499 idxt = num_entries + 1
1500 do idx=1,num_entries
1501 if ((index(fieldname(idx), trim(short_fieldnm)) /= 0) .and. &
1502 (len_trim(fieldname(idx)) == len_trim(short_fieldnm))) then
1503 if (index(input_name,trim(from_input(idx))) /= 0 .or. &
1504 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1505 idxt = idx
1506 end if
1507 end if
1508 end do
1509 idx = idxt
1510 if (idx > num_entries) then
1511 call mprintf(.true.,WARN,'Entry in METGRID.TBL not found for field %s. '// &
1512 'Default options will be used for this field!', s1=short_fieldnm)
1513 idx = num_entries ! The last entry is a default
1514 end if
1516 if (process_bdy_width == 0) then
1517 process_width = max(ed1-sd1+1, ed2-sd2+1)
1518 else
1519 process_width = process_bdy_width
1520 ! Add two extra rows/cols to accommodate staggered fields: one extra row/col for
1521 ! averaging to mass points in real, and one beyond that for averaging during
1522 ! wind rotation
1523 if (ifieldstagger /= M) process_width = process_width + 2
1524 end if
1526 field%header%dim1(1) = sm1
1527 field%header%dim1(2) = em1
1528 field%header%dim2(1) = sm2
1529 field%header%dim2(2) = em2
1530 field%map%stagger = ifieldstagger
1531 if (.not. associated(field%r_arr)) then
1532 allocate(field%r_arr(sm1:em1,sm2:em2))
1533 end if
1535 interp_mask_status = 1
1536 interp_land_mask_status = 1
1537 interp_water_mask_status = 1
1539 if (interp_mask(idx) /= ' ') then
1540 mask_field%header%version = 1
1541 mask_field%header%forecast_hour = 0.
1542 mask_field%header%field = trim(interp_mask(idx))//'.mask'
1543 mask_field%header%vertical_coord = 'none'
1544 mask_field%header%vertical_level = 1
1546 call storage_get_field(mask_field, interp_mask_status)
1548 end if
1549 if (interp_land_mask(idx) /= ' ') then
1550 mask_land_field%header%version = 1
1551 mask_land_field%header%forecast_hour = 0.
1552 mask_land_field%header%field = trim(interp_land_mask(idx))//'.mask'
1553 mask_land_field%header%vertical_coord = 'none'
1554 mask_land_field%header%vertical_level = 1
1556 call storage_get_field(mask_land_field, interp_land_mask_status)
1558 end if
1559 if (interp_water_mask(idx) /= ' ') then
1560 mask_water_field%header%version = 1
1561 mask_water_field%header%forecast_hour = 0.
1562 mask_water_field%header%field = trim(interp_water_mask(idx))//'.mask'
1563 mask_water_field%header%vertical_coord = 'none'
1564 mask_water_field%header%vertical_level = 1
1566 call storage_get_field(mask_water_field, interp_water_mask_status)
1568 end if
1570 interp_array => interp_array_from_string(interp_method(idx))
1571 interp_opts => interp_options_from_string(interp_method(idx))
1574 !
1575 ! Interpolate using average_gcell interpolation method
1576 !
1577 if (do_gcell_interp) then
1578 !BPR BEGIN
1579 !If a lower priority source of the current field has already been read
1580 !in, the results of that input are currently in field%r_arr
1581 !Pass COPY of field%r_arr into accum_continous because in accum_continuous
1582 !it will set the input variable to zero over points covered by the
1583 !current source and then determine the appropriate value to place at
1584 !that point. This will overwrite data already put in field%r_arr by
1585 !lower priority sources.
1586 !This is problematic if the current source results in missing values
1587 !over parts of the area covered by the current source where a lower
1588 !priority source has already provided a value. In this case, if one
1589 !passes in field%r_arr, it will overwrite the value provided by the
1590 !lower priority source with zero.
1591 r_arr_cur_source = field%r_arr
1592 !BPR END
1593 allocate(data_count(sm1:em1,sm2:em2))
1594 data_count = 0.
1596 if (interp_mask_status == 0) then
1597 call accum_continuous(slab, &
1598 minx, maxx, miny, maxy, 1, 1, bdr, &
1599 !BPR BEGIN
1600 !Pass COPY of field%r_arr instead of field%r_arr itself
1601 !field%r_arr, data_count, &
1602 r_arr_cur_source, data_count, &
1603 !BPR END
1604 sm1, em1, sm2, em2, 1, 1, &
1605 istagger, &
1606 new_pts, missing_value(idx), interp_mask_val(idx), interp_mask_relational(idx), mask_field%r_arr)
1607 else
1608 call accum_continuous(slab, &
1609 minx, maxx, miny, maxy, 1, 1, bdr, &
1610 !BPR BEGIN
1611 !Pass COPY of field%r_arr instead of field%r_arr itself
1612 !field%r_arr, data_count, &
1613 r_arr_cur_source, data_count, &
1614 !BPR END
1615 sm1, em1, sm2, em2, 1, 1, &
1616 istagger, &
1617 new_pts, missing_value(idx), -1.) ! The -1 is the maskval, but since we
1618 ! we do not give an optional mask, no
1619 ! no need to worry about -1s in data
1620 end if
1622 orig_selected_proj = iget_selected_domain()
1623 call select_domain(SOURCE_PROJ)
1624 do j=sm2,em2
1625 do i=sm1,em1
1627 if (abs(i - sd1) >= process_width .and. (abs(i - ed1) >= process_width) .and. &
1628 abs(j - sd2) >= process_width .and. (abs(j - ed2) >= process_width)) then
1629 field%r_arr(i,j) = fill_missing(idx)
1630 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1631 cycle
1632 end if
1634 if (present(landmask)) then
1636 if (landmask(i,j) /= masked(idx)) then
1637 if (data_count(i,j) > 0.) then
1638 !BPR BEGIN
1639 !accum_continuous is now passed a copy of field%r_arr (r_arr_cur_source)
1640 !instead of field%r_arr itself so that it does not set
1641 !field%r_arr to zero where the input source is marked as missing
1642 !field%r_arr(i,j) = field%r_arr(i,j) / data_count(i,j)
1643 field%r_arr(i,j) = r_arr_cur_source(i,j) / data_count(i,j)
1644 !BPR END
1645 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1646 else
1648 if (interp_mask_status == 0) then
1649 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1650 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1651 mask_relational=interp_mask_relational(idx), &
1652 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1653 else
1654 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1655 minx, maxx, miny, maxy, bdr, missing_value(idx))
1656 end if
1658 if (temp /= missing_value(idx)) then
1659 field%r_arr(i,j) = temp
1660 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1661 end if
1663 end if
1664 else
1665 field%r_arr(i,j) = fill_missing(idx)
1666 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1667 end if
1669 if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
1670 .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
1671 field%r_arr(i,j) = fill_missing(idx)
1673 ! Assume that if missing fill value is other than default, then user has asked
1674 ! to fill in any missing values, and we can consider this point to have
1675 ! received a valid value
1676 if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1677 end if
1679 else
1681 if (data_count(i,j) > 0.) then
1682 !BPR BEGIN
1683 !accum_continuous is now passed a copy of field%r_arr (r_arr_cur_source)
1684 !instead of field%r_arr itself so that it does not set
1685 !field%r_arr to zero where the input source is marked as missing
1686 !field%r_arr(i,j) = field%r_arr(i,j) / data_count(i,j)
1687 field%r_arr(i,j) = r_arr_cur_source(i,j) / data_count(i,j)
1688 !BPR END
1689 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1690 else
1692 if (interp_mask_status == 0) then
1693 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1694 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1695 mask_relational=interp_mask_relational(idx), &
1696 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1697 else
1698 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1699 minx, maxx, miny, maxy, bdr, missing_value(idx))
1700 end if
1702 if (temp /= missing_value(idx)) then
1703 field%r_arr(i,j) = temp
1704 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1705 end if
1707 end if
1709 if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
1710 .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
1711 field%r_arr(i,j) = fill_missing(idx)
1713 ! Assume that if missing fill value is other than default, then user has asked
1714 ! to fill in any missing values, and we can consider this point to have
1715 ! received a valid value
1716 if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1717 end if
1719 end if
1721 end do
1722 end do
1723 call select_domain(orig_selected_proj)
1724 deallocate(data_count)
1726 !
1727 ! No average_gcell interpolation method
1728 !
1729 else
1731 orig_selected_proj = iget_selected_domain()
1732 call select_domain(SOURCE_PROJ)
1733 do j=sm2,em2
1734 do i=sm1,em1
1736 if (abs(i - sd1) >= process_width .and. (abs(i - ed1) >= process_width) .and. &
1737 abs(j - sd2) >= process_width .and. (abs(j - ed2) >= process_width)) then
1738 field%r_arr(i,j) = fill_missing(idx)
1739 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1740 cycle
1741 end if
1743 if (present(landmask)) then
1745 if (masked(idx) == MASKED_BOTH) then
1747 if (landmask(i,j) == 0) then ! WATER POINT
1749 if (interp_land_mask_status == 0) then
1750 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1751 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1752 mask_relational=interp_land_mask_relational(idx), &
1753 mask_val=interp_land_mask_val(idx), mask_field=mask_land_field%r_arr)
1754 else
1755 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1756 minx, maxx, miny, maxy, bdr, missing_value(idx))
1757 end if
1759 else if (landmask(i,j) == 1) then ! LAND POINT
1761 if (interp_water_mask_status == 0) then
1762 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1763 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1764 mask_relational=interp_water_mask_relational(idx), &
1765 mask_val=interp_water_mask_val(idx), mask_field=mask_water_field%r_arr)
1766 else
1767 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1768 minx, maxx, miny, maxy, bdr, missing_value(idx))
1769 end if
1771 end if
1773 else if (landmask(i,j) /= masked(idx)) then
1775 if (interp_mask_status == 0) then
1776 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1777 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1778 mask_relational=interp_mask_relational(idx), &
1779 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1780 else
1781 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1782 minx, maxx, miny, maxy, bdr, missing_value(idx))
1783 end if
1785 else
1786 temp = missing_value(idx)
1787 end if
1789 ! No landmask for this field
1790 else
1792 if (interp_mask_status == 0) then
1793 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1794 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1795 mask_relational=interp_mask_relational(idx), &
1796 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1797 else
1798 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1799 minx, maxx, miny, maxy, bdr, missing_value(idx))
1800 end if
1802 end if
1804 if (temp /= missing_value(idx)) then
1805 field%r_arr(i,j) = temp
1806 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1807 else if (present(landmask)) then
1808 if (landmask(i,j) == masked(idx)) then
1809 field%r_arr(i,j) = fill_missing(idx)
1810 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1811 end if
1812 end if
1814 if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
1815 .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
1816 field%r_arr(i,j) = fill_missing(idx)
1818 ! Assume that if missing fill value is other than default, then user has asked
1819 ! to fill in any missing values, and we can consider this point to have
1820 ! received a valid value
1821 if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1822 end if
1824 end do
1825 end do
1826 call select_domain(orig_selected_proj)
1827 end if
1829 deallocate(interp_array)
1830 deallocate(interp_opts)
1832 end subroutine interp_met_field
1835 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1836 ! Name: interp_to_latlon
1837 !
1838 ! Purpose:
1839 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1840 function interp_to_latlon(rlat, rlon, istagger, interp_method_list, interp_opt_list, slab, &
1841 minx, maxx, miny, maxy, bdr, source_missing_value, &
1842 mask_field, mask_relational, mask_val)
1844 use interp_module
1845 use llxy_module
1847 implicit none
1849 ! Arguments
1850 integer, intent(in) :: minx, maxx, miny, maxy, bdr, istagger
1851 integer, dimension(:), intent(in) :: interp_method_list
1852 integer, dimension(:), intent(in) :: interp_opt_list
1853 real, intent(in) :: rlat, rlon, source_missing_value
1854 real, dimension(minx:maxx,miny:maxy), intent(in) :: slab
1855 real, intent(in), optional :: mask_val
1856 real, dimension(minx:maxx,miny:maxy), intent(in), optional :: mask_field
1857 character(len=1), intent(in), optional :: mask_relational
1859 ! Return value
1860 real :: interp_to_latlon
1862 ! Local variables
1863 real :: rx, ry
1865 interp_to_latlon = source_missing_value
1867 call lltoxy(rlat, rlon, rx, ry, istagger)
1868 if (rx >= minx+bdr-0.5 .and. rx <= maxx-bdr+0.5) then
1869 if (present(mask_field) .and. present(mask_val) .and. present (mask_relational)) then
1870 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
1871 interp_method_list, interp_opt_list, 1, mask_relational, mask_val, mask_field)
1872 else if (present(mask_field) .and. present(mask_val)) then
1873 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
1874 interp_method_list, interp_opt_list, 1, maskval=mask_val, mask_array=mask_field)
1875 else
1876 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
1877 interp_method_list, interp_opt_list, 1)
1878 end if
1879 else
1880 interp_to_latlon = source_missing_value
1881 end if
1883 if (interp_to_latlon == source_missing_value) then
1885 ! Try a lon in the range 0. to 360.; all lons in the xlon
1886 ! array should be in the range -180. to 180.
1887 if (rlon < 0.) then
1888 call lltoxy(rlat, rlon+360., rx, ry, istagger)
1889 if (rx >= minx+bdr-0.5 .and. rx <= maxx-bdr+0.5) then
1890 if (present(mask_field) .and. present(mask_val) .and. present(mask_relational)) then
1891 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
1892 1, 1, source_missing_value, &
1893 interp_method_list, interp_opt_list, 1, &
1894 mask_relational, mask_val, mask_field)
1895 else if (present(mask_field) .and. present(mask_val)) then
1896 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
1897 1, 1, source_missing_value, &
1898 interp_method_list, interp_opt_list, 1, &
1899 maskval=mask_val, mask_array=mask_field)
1900 else
1901 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
1902 1, 1, source_missing_value, &
1903 interp_method_list, interp_opt_list, 1)
1904 end if
1905 else
1906 interp_to_latlon = source_missing_value
1907 end if
1909 end if
1911 end if
1913 return
1915 end function interp_to_latlon
1918 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1919 ! Name: get_bottom_top_dim
1920 !
1921 ! Purpose:
1922 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1923 subroutine get_bottom_top_dim(bottom_top_dim)
1925 use interp_option_module
1926 use list_module
1927 use storage_module
1929 implicit none
1931 ! Arguments
1932 integer, intent(out) :: bottom_top_dim
1934 ! Local variables
1935 integer :: i, j
1936 integer, pointer, dimension(:) :: field_levels
1937 character (len=32) :: z_dim
1938 type (fg_input), pointer, dimension(:) :: headers
1939 type (list) :: temp_levels
1941 !CWH Initialize local pointer variables
1942 nullify(field_levels)
1943 nullify(headers)
1945 ! Initialize a list to store levels that are found for 3-d fields
1946 call list_init(temp_levels)
1948 ! Get a list of all time-dependent fields (given by their headers) from
1949 ! the storage module
1950 call storage_get_td_headers(headers)
1952 !
1953 ! Given headers of all fields, we first build a list of all possible levels
1954 ! for 3-d met fields (excluding sea-level, though).
1955 !
1956 do i=1,size(headers)
1957 call get_z_dim_name(headers(i)%header%field, z_dim)
1959 ! We only want to consider 3-d met fields
1960 if (z_dim(1:18) == 'num_metgrid_levels') then
1962 ! Find out what levels the current field has
1963 call storage_get_levels(headers(i), field_levels)
1964 do j=1,size(field_levels)
1966 ! If this level has not yet been encountered, add it to our list
1967 if (.not. list_search(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))) then
1968 if (field_levels(j) /= 201300) then
1969 call list_insert(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))
1970 end if
1971 end if
1973 end do
1975 deallocate(field_levels)
1977 end if
1979 end do
1981 bottom_top_dim = list_length(temp_levels)
1983 call list_destroy(temp_levels)
1984 deallocate(headers)
1986 end subroutine get_bottom_top_dim
1989 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1990 ! Name: fill_missing_levels
1991 !
1992 ! Purpose:
1993 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1994 subroutine fill_missing_levels(output_flags)
1996 use interp_option_module
1997 use list_module
1998 use module_debug
1999 use module_mergesort
2000 use storage_module
2002 implicit none
2004 ! Arguments
2005 character (len=128), dimension(:), intent(inout) :: output_flags
2007 ! Local variables
2008 integer :: i, ii, j, ix, jx, k, lower, upper, temp, istatus
2009 integer, pointer, dimension(:) :: union_levels, field_levels
2010 real, pointer, dimension(:) :: r_union_levels
2011 character (len=128) :: clevel
2012 type (fg_input) :: lower_field, upper_field, new_field, search_field
2013 type (fg_input), pointer, dimension(:) :: headers, all_headers
2014 type (list) :: temp_levels
2015 type (list_item), pointer, dimension(:) :: keys
2017 ! CWH Initialize local pointer variables
2018 nullify(union_levels)
2019 nullify(field_levels)
2020 nullify(r_union_levels)
2021 nullify(headers)
2022 nullify(all_headers)
2023 nullify(keys)
2025 ! Initialize a list to store levels that are found for 3-d fields
2026 call list_init(temp_levels)
2028 ! Get a list of all fields (given by their headers) from the storage module
2029 call storage_get_td_headers(headers)
2030 call storage_get_all_headers(all_headers)
2032 !
2033 ! Given headers of all fields, we first build a list of all possible levels
2034 ! for 3-d met fields (excluding sea-level, though).
2035 !
2036 do i=1,size(headers)
2038 ! Find out what levels the current field has
2039 call storage_get_levels(headers(i), field_levels)
2040 do j=1,size(field_levels)
2042 ! If this level has not yet been encountered, add it to our list
2043 if (.not. list_search(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))) then
2044 if (field_levels(j) /= 201300) then
2045 call list_insert(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))
2046 end if
2047 end if
2049 end do
2051 deallocate(field_levels)
2053 end do
2055 if (list_length(temp_levels) > 0) then
2057 !
2058 ! With all possible levels stored in a list, get an array of levels, sorted
2059 ! in decreasing order
2060 !
2061 i = 0
2062 allocate(union_levels(list_length(temp_levels)))
2063 do while (list_length(temp_levels) > 0)
2064 i = i + 1
2065 call list_get_first_item(temp_levels, ikey=union_levels(i), ivalue=temp)
2066 end do
2067 call mergesort(union_levels, 1, size(union_levels))
2069 allocate(r_union_levels(size(union_levels)))
2070 do i=1,size(union_levels)
2071 r_union_levels(i) = real(union_levels(i))
2072 end do
2074 !
2075 ! With a sorted, complete list of levels, we need
2076 ! to go back and fill in missing levels for each 3-d field
2077 !
2078 do i=1,size(headers)
2080 !
2081 ! Find entry in METGRID.TBL for this field, if one exists; if it does, then the
2082 ! entry may tell us how to get values for the current field at the missing level
2083 !
2084 do ii=1,num_entries
2085 if (fieldname(ii) == headers(i)%header%field) exit
2086 end do
2087 if (ii <= num_entries) then
2088 call dup(headers(i),new_field)
2089 nullify(new_field%valid_mask)
2090 nullify(new_field%modified_mask)
2091 call fill_field(new_field, ii, output_flags, r_union_levels)
2092 end if
2094 end do
2096 deallocate(union_levels)
2097 deallocate(r_union_levels)
2098 deallocate(headers)
2100 call storage_get_td_headers(headers)
2102 !
2103 ! Now we may need to vertically interpolate to missing values in 3-d fields
2104 !
2105 do i=1,size(headers)
2107 call storage_get_levels(headers(i), field_levels)
2109 ! If this isn't a 3-d array, nothing to do
2110 if (size(field_levels) > 1) then
2112 do k=1,size(field_levels)
2113 call dup(headers(i),search_field)
2114 search_field%header%vertical_level = field_levels(k)
2115 call storage_get_field(search_field,istatus)
2116 if (istatus == 0) then
2117 JLOOP: do jx=search_field%header%dim2(1),search_field%header%dim2(2)
2118 ILOOP: do ix=search_field%header%dim1(1),search_field%header%dim1(2)
2119 if (.not. bitarray_test(search_field%valid_mask, &
2120 ix-search_field%header%dim1(1)+1, &
2121 jx-search_field%header%dim2(1)+1)) then
2123 call dup(search_field, lower_field)
2124 do lower=k-1,1,-1
2125 lower_field%header%vertical_level = field_levels(lower)
2126 call storage_get_field(lower_field,istatus)
2127 if (bitarray_test(lower_field%valid_mask, &
2128 ix-search_field%header%dim1(1)+1, &
2129 jx-search_field%header%dim2(1)+1)) &
2130 exit
2132 end do
2134 call dup(search_field, upper_field)
2135 do upper=k+1,size(field_levels)
2136 upper_field%header%vertical_level = field_levels(upper)
2137 call storage_get_field(upper_field,istatus)
2138 if (bitarray_test(upper_field%valid_mask, &
2139 ix-search_field%header%dim1(1)+1, &
2140 jx-search_field%header%dim2(1)+1)) &
2141 exit
2143 end do
2144 if (upper <= size(field_levels) .and. lower >= 1) then
2145 search_field%r_arr(ix,jx) = real(abs(field_levels(upper)-field_levels(k))) &
2146 / real(abs(field_levels(upper)-field_levels(lower))) &
2147 * lower_field%r_arr(ix,jx) &
2148 + real(abs(field_levels(k)-field_levels(lower))) &
2149 / real(abs(field_levels(upper)-field_levels(lower))) &
2150 * upper_field%r_arr(ix,jx)
2151 call bitarray_set(search_field%valid_mask, &
2152 ix-search_field%header%dim1(1)+1, &
2153 jx-search_field%header%dim2(1)+1)
2154 end if
2155 end if
2156 end do ILOOP
2157 end do JLOOP
2158 else
2159 call mprintf(.true.,ERROR, &
2160 'This is bad, could not get %s at level %i.', &
2161 s1=trim(search_field%header%field), i1=field_levels(k))
2162 end if
2163 end do
2165 end if
2167 deallocate(field_levels)
2169 end do
2171 end if
2173 call list_destroy(temp_levels)
2174 deallocate(all_headers)
2175 deallocate(headers)
2177 end subroutine fill_missing_levels
2180 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2181 ! Name: create_derived_fields
2182 !
2183 ! Purpose:
2184 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2185 subroutine create_derived_fields(arg_gridtype, hdate, xfcst, &
2186 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
2187 we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e, &
2188 created_this_field, output_flags)
2190 use interp_option_module
2191 use list_module
2192 use module_mergesort
2193 use storage_module
2195 implicit none
2197 ! Arguments
2198 integer, intent(in) :: we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
2199 we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e
2200 real, intent(in) :: xfcst
2201 logical, dimension(:), intent(inout) :: created_this_field
2202 character (len=1), intent(in) :: arg_gridtype
2203 character (len=24), intent(in) :: hdate
2204 character (len=128), dimension(:), intent(inout) :: output_flags
2206 ! Local variables
2207 integer :: idx, i, j, istatus
2208 type (fg_input) :: field
2210 ! Initialize fg_input structure to store the field
2211 field%header%version = 1
2212 field%header%date = hdate//' '
2213 field%header%time_dependent = .true.
2214 field%header%mask_field = .false.
2215 field%header%constant_field = .false.
2216 field%header%forecast_hour = xfcst
2217 field%header%fg_source = 'Derived from FG'
2218 field%header%field = ' '
2219 field%header%units = ' '
2220 field%header%description = ' '
2221 field%header%vertical_level = 0
2222 field%header%sr_x = 1
2223 field%header%sr_y = 1
2224 field%header%array_order = 'XY '
2225 field%header%is_wind_grid_rel = .true.
2226 field%header%array_has_missing_values = .false.
2227 nullify(field%r_arr)
2228 nullify(field%valid_mask)
2229 nullify(field%modified_mask)
2231 !
2232 ! Check each entry in METGRID.TBL to see whether it is a derive field
2233 !
2234 do idx=1,num_entries
2235 if (is_derived_field(idx)) then
2237 created_this_field(idx) = .true.
2239 call mprintf(.true.,INFORM,'Going to create the field %s',s1=fieldname(idx))
2241 ! Intialize more fields in storage structure
2242 field%header%field = fieldname(idx)
2243 call get_z_dim_name(fieldname(idx),field%header%vertical_coord)
2244 field%map%stagger = output_stagger(idx)
2245 if (arg_gridtype == 'E') then
2246 field%header%dim1(1) = we_mem_s
2247 field%header%dim1(2) = we_mem_e
2248 field%header%dim2(1) = sn_mem_s
2249 field%header%dim2(2) = sn_mem_e
2250 else if (arg_gridtype == 'C') then
2251 if (output_stagger(idx) == M) then
2252 field%header%dim1(1) = we_mem_s
2253 field%header%dim1(2) = we_mem_e
2254 field%header%dim2(1) = sn_mem_s
2255 field%header%dim2(2) = sn_mem_e
2256 else if (output_stagger(idx) == U) then
2257 field%header%dim1(1) = we_mem_stag_s
2258 field%header%dim1(2) = we_mem_stag_e
2259 field%header%dim2(1) = sn_mem_s
2260 field%header%dim2(2) = sn_mem_e
2261 else if (output_stagger(idx) == V) then
2262 field%header%dim1(1) = we_mem_s
2263 field%header%dim1(2) = we_mem_e
2264 field%header%dim2(1) = sn_mem_stag_s
2265 field%header%dim2(2) = sn_mem_stag_e
2266 end if
2267 end if
2269 call fill_field(field, idx, output_flags)
2271 end if
2272 end do
2275 end subroutine create_derived_fields
2278 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2279 ! Name: fill_field
2280 !
2281 ! Purpose:
2282 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2283 subroutine fill_field(field, idx, output_flags, all_level_list)
2285 use interp_option_module
2286 use list_module
2287 use module_mergesort
2288 use storage_module
2290 implicit none
2292 ! Arguments
2293 integer, intent(in) :: idx
2294 type (fg_input), intent(inout) :: field
2295 character (len=128), dimension(:), intent(inout) :: output_flags
2296 real, dimension(:), intent(in), optional :: all_level_list
2298 ! Local variables
2299 integer :: i, j, istatus, isrclevel
2300 integer, pointer, dimension(:) :: all_list
2301 real :: rfillconst, rlevel, rsrclevel
2302 type (fg_input) :: query_field
2303 type (list_item), pointer, dimension(:) :: keys
2304 character (len=128) :: asrcname
2305 logical :: filled_all_lev
2307 !CWH Initialize local pointer variables
2308 nullify(all_list)
2309 nullify(keys)
2311 filled_all_lev = .false.
2313 !
2314 ! Get a list of all levels to be filled for this field
2315 !
2316 keys => list_get_keys(fill_lev_list(idx))
2318 do i=1,list_length(fill_lev_list(idx))
2320 !
2321 ! First handle a specification for levels "all"
2322 !
2323 if (trim(keys(i)%ckey) == 'all') then
2325 ! We only want to fill all levels if we haven't already filled "all" of them
2326 if (.not. filled_all_lev) then
2328 filled_all_lev = .true.
2330 query_field%header%time_dependent = .true.
2331 query_field%header%field = ' '
2332 nullify(query_field%r_arr)
2333 nullify(query_field%valid_mask)
2334 nullify(query_field%modified_mask)
2336 ! See if we are filling this level with a constant
2337 call get_constant_fill_lev(keys(i)%cvalue, rfillconst, istatus)
2338 if (istatus == 0) then
2339 if (present(all_level_list)) then
2340 do j=1,size(all_level_list)
2341 call create_level(field, real(all_level_list(j)), idx, output_flags, rfillconst=rfillconst)
2342 end do
2343 else
2344 query_field%header%field = level_template(idx)
2345 nullify(all_list)
2346 call storage_get_levels(query_field, all_list)
2347 if (associated(all_list)) then
2348 do j=1,size(all_list)
2349 call create_level(field, real(all_list(j)), idx, output_flags, rfillconst=rfillconst)
2350 end do
2351 deallocate(all_list)
2352 end if
2353 end if
2355 ! Else see if we are filling this level with a constant equal
2356 ! to the value of the level
2357 else if (trim(keys(i)%cvalue) == 'vertical_index') then
2358 if (present(all_level_list)) then
2359 do j=1,size(all_level_list)
2360 call create_level(field, real(all_level_list(j)), idx, output_flags, &
2361 rfillconst=real(all_level_list(j)))
2362 end do
2363 else
2364 query_field%header%field = level_template(idx)
2365 nullify(all_list)
2366 call storage_get_levels(query_field, all_list)
2367 if (associated(all_list)) then
2368 do j=1,size(all_list)
2369 call create_level(field, real(all_list(j)), idx, output_flags, rfillconst=real(all_list(j)))
2370 end do
2371 deallocate(all_list)
2372 end if
2373 end if
2375 ! Else, we assume that it is a field from which we are copying levels
2376 else
2377 if (present(all_level_list)) then
2378 do j=1,size(all_level_list)
2379 call create_level(field, real(all_level_list(j)), idx, output_flags, &
2380 asrcname=keys(i)%cvalue, rsrclevel=real(all_level_list(j)))
2381 end do
2382 else
2383 query_field%header%field = keys(i)%cvalue ! Use same levels as source field, not level_template
2384 nullify(all_list)
2385 call storage_get_levels(query_field, all_list)
2386 if (associated(all_list)) then
2387 do j=1,size(all_list)
2388 call create_level(field, real(all_list(j)), idx, output_flags, &
2389 asrcname=keys(i)%cvalue, rsrclevel=real(all_list(j)))
2390 end do
2391 deallocate(all_list)
2393 else
2395 ! If the field doesn't have any levels (or does not exist) then we have not
2396 ! really filled all levels at this point.
2397 filled_all_lev = .false.
2398 end if
2399 end if
2401 end if
2402 end if
2404 !
2405 ! Handle individually specified levels
2406 !
2407 else
2409 read(keys(i)%ckey,*) rlevel
2411 ! See if we are filling this level with a constant
2412 call get_constant_fill_lev(keys(i)%cvalue, rfillconst, istatus)
2413 if (istatus == 0) then
2414 call create_level(field, rlevel, idx, output_flags, rfillconst=rfillconst)
2416 ! Otherwise, we are filling from another level
2417 else
2418 call get_fill_src_level(keys(i)%cvalue, asrcname, isrclevel)
2419 rsrclevel = real(isrclevel)
2420 call create_level(field, rlevel, idx, output_flags, &
2421 asrcname=asrcname, rsrclevel=rsrclevel)
2423 end if
2424 end if
2425 end do
2427 if (associated(keys)) deallocate(keys)
2429 end subroutine fill_field
2432 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2433 ! Name: create_level
2434 !
2435 ! Purpose:
2436 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2437 subroutine create_level(field_template, rlevel, idx, output_flags, &
2438 rfillconst, asrcname, rsrclevel)
2440 use storage_module
2441 use interp_option_module
2443 implicit none
2445 ! Arguments
2446 type (fg_input), intent(inout) :: field_template
2447 real, intent(in) :: rlevel
2448 integer, intent(in) :: idx
2449 character (len=128), dimension(:), intent(inout) :: output_flags
2450 real, intent(in), optional :: rfillconst, rsrclevel
2451 character (len=128), intent(in), optional :: asrcname
2453 ! Local variables
2454 integer :: i, j, istatus
2455 integer :: sm1,em1,sm2,em2
2456 type (fg_input) :: query_field
2458 !
2459 ! Check to make sure optional arguments are sane
2460 !
2461 if (present(rfillconst) .and. (present(asrcname) .or. present(rsrclevel))) then
2462 call mprintf(.true.,ERROR,'A call to create_level() cannot be given specifications '// &
2463 'for both a constant fill value and a source level.')
2465 else if ((present(asrcname) .and. .not. present(rsrclevel)) .or. &
2466 (.not. present(asrcname) .and. present(rsrclevel))) then
2467 call mprintf(.true.,ERROR,'Neither or both of optional arguments asrcname and '// &
2468 'rsrclevel must be specified to subroutine create_level().')
2470 else if (.not. present(rfillconst) .and. &
2471 .not. present(asrcname) .and. &
2472 .not. present(rsrclevel)) then
2473 call mprintf(.true.,ERROR,'A call to create_level() must be given either a specification '// &
2474 'for a constant fill value or a source level.')
2475 end if
2477 query_field%header%time_dependent = .true.
2478 query_field%header%field = field_template%header%field
2479 query_field%header%vertical_level = rlevel
2480 nullify(query_field%r_arr)
2481 nullify(query_field%valid_mask)
2482 nullify(query_field%modified_mask)
2484 call storage_query_field(query_field, istatus)
2485 if (istatus == 0) then
2486 call mprintf(.true.,INFORM,'%s at level %f already exists; leaving it alone.', &
2487 s1=field_template%header%field, f1=rlevel)
2488 return
2489 end if
2491 sm1 = field_template%header%dim1(1)
2492 em1 = field_template%header%dim1(2)
2493 sm2 = field_template%header%dim2(1)
2494 em2 = field_template%header%dim2(2)
2496 !
2497 ! Handle constant fill value case
2498 !
2499 if (present(rfillconst)) then
2501 field_template%header%vertical_level = rlevel
2502 allocate(field_template%r_arr(sm1:em1,sm2:em2))
2503 allocate(field_template%valid_mask)
2504 allocate(field_template%modified_mask)
2505 call bitarray_create(field_template%valid_mask, em1-sm1+1, em2-sm2+1)
2506 call bitarray_create(field_template%modified_mask, em1-sm1+1, em2-sm2+1)
2508 field_template%r_arr = rfillconst
2510 do j=sm2,em2
2511 do i=sm1,em1
2512 call bitarray_set(field_template%valid_mask, i-sm1+1, j-sm2+1)
2513 end do
2514 end do
2516 call storage_put_field(field_template)
2518 if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
2519 output_flags(idx) = flag_in_output(idx)
2520 end if
2522 !
2523 ! Handle source field and source level case
2524 !
2525 else if (present(asrcname) .and. present(rsrclevel)) then
2527 query_field%header%field = ' '
2528 query_field%header%field = asrcname
2529 query_field%header%vertical_level = rsrclevel
2531 ! Check to see whether the requested source field exists at the requested level
2532 call storage_query_field(query_field, istatus)
2534 if (istatus == 0) then
2536 ! Read in requested field at requested level
2537 call storage_get_field(query_field, istatus)
2538 if ((query_field%header%dim1(1) /= field_template%header%dim1(1)) .or. &
2539 (query_field%header%dim1(2) /= field_template%header%dim1(2)) .or. &
2540 (query_field%header%dim2(1) /= field_template%header%dim2(1)) .or. &
2541 (query_field%header%dim2(2) /= field_template%header%dim2(2))) then
2542 call mprintf(.true.,ERROR,'Dimensions for %s do not match those of %s. This is '// &
2543 'probably because the staggerings of the fields do not match.', &
2544 s1=query_field%header%field, s2=field_template%header%field)
2545 end if
2547 field_template%header%vertical_level = rlevel
2548 allocate(field_template%r_arr(sm1:em1,sm2:em2))
2549 allocate(field_template%valid_mask)
2550 allocate(field_template%modified_mask)
2551 call bitarray_create(field_template%valid_mask, em1-sm1+1, em2-sm2+1)
2552 call bitarray_create(field_template%modified_mask, em1-sm1+1, em2-sm2+1)
2554 field_template%r_arr = query_field%r_arr
2556 ! We should retain information about which points in the field are valid
2557 do j=sm2,em2
2558 do i=sm1,em1
2559 if (bitarray_test(query_field%valid_mask, i-sm1+1, j-sm2+1)) then
2560 call bitarray_set(field_template%valid_mask, i-sm1+1, j-sm2+1)
2561 end if
2562 end do
2563 end do
2565 call storage_put_field(field_template)
2567 if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
2568 output_flags(idx) = flag_in_output(idx)
2569 end if
2571 else
2572 call mprintf(.true.,INFORM,'Couldn''t find %s at level %f to fill level %f of %s.', &
2573 s1=asrcname,f1=rsrclevel,f2=rlevel,s2=field_template%header%field)
2574 end if
2576 end if
2578 end subroutine create_level
2581 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2582 ! Name: accum_continuous
2583 !
2584 ! Purpose: Sum up all of the source data points whose nearest neighbor in the
2585 ! model grid is the specified model grid point.
2586 !
2587 ! NOTE: When processing the source tile, those source points that are
2588 ! closest to a different model grid point will be added to the totals for
2589 ! such grid points; thus, an entire source tile will be processed at a time.
2590 ! This routine really processes for all model grid points that are
2591 ! within a source tile, and not just for a single grid point.
2592 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2593 subroutine accum_continuous(src_array, &
2594 src_min_x, src_max_x, src_min_y, src_max_y, src_min_z, src_max_z, bdr_width, &
2595 dst_array, n, &
2596 start_i, end_i, start_j, end_j, start_k, end_k, &
2597 istagger, &
2598 new_pts, msgval, maskval, mask_relational, mask_array, sr_x, sr_y)
2600 use bitarray_module
2601 use misc_definitions_module
2603 implicit none
2605 ! Arguments
2606 integer, intent(in) :: start_i, end_i, start_j, end_j, start_k, end_k, istagger, &
2607 src_min_x, src_max_x, src_min_y, src_max_y, src_min_z, src_max_z, bdr_width
2608 real, intent(in) :: maskval, msgval
2609 real, dimension(src_min_x:src_max_x, src_min_y:src_max_y, src_min_z:src_max_z), intent(in) :: src_array
2610 real, dimension(start_i:end_i, start_j:end_j, start_k:end_k), intent(inout) :: dst_array, n
2611 real, dimension(src_min_x:src_max_x, src_min_y:src_max_y), intent(in), optional :: mask_array
2612 integer, intent(in), optional :: sr_x, sr_y
2613 type (bitarray), intent(inout) :: new_pts
2614 character(len=1), intent(in), optional :: mask_relational
2616 ! Local variables
2617 integer :: i, j
2618 integer, pointer, dimension(:,:,:) :: where_maps_to
2619 real :: rsr_x, rsr_y
2621 rsr_x = 1.0
2622 rsr_y = 1.0
2623 if (present(sr_x)) rsr_x = real(sr_x)
2624 if (present(sr_y)) rsr_y = real(sr_y)
2626 allocate(where_maps_to(src_min_x:src_max_x,src_min_y:src_max_y,2))
2627 do i=src_min_x,src_max_x
2628 do j=src_min_y,src_max_y
2629 where_maps_to(i,j,1) = NOT_PROCESSED
2630 end do
2631 end do
2633 call process_continuous_block(src_array, where_maps_to, &
2634 src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
2635 src_min_x+bdr_width, src_min_y, src_min_z, &
2636 src_max_x-bdr_width, src_max_y, src_max_z, &
2637 dst_array, n, start_i, end_i, start_j, end_j, start_k, end_k, &
2638 istagger, &
2639 new_pts, rsr_x, rsr_y, msgval, maskval, mask_relational, mask_array)
2641 deallocate(where_maps_to)
2643 end subroutine accum_continuous
2646 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2647 ! Name: process_continuous_block
2648 !
2649 ! Purpose: To recursively process a subarray of continuous data, adding the
2650 ! points in a block to the sum for their nearest grid point. The nearest
2651 ! neighbor may be estimated in some cases; for example, if the four corners
2652 ! of a subarray all have the same nearest grid point, all elements in the
2653 ! subarray are added to that grid point.
2654 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2655 recursive subroutine process_continuous_block(tile_array, where_maps_to, &
2656 src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
2657 min_i, min_j, min_k, max_i, max_j, max_k, &
2658 dst_array, n, &
2659 start_x, end_x, start_y, end_y, start_z, end_z, &
2660 istagger, &
2661 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2663 use bitarray_module
2664 use llxy_module
2665 use misc_definitions_module
2667 implicit none
2669 ! Arguments
2670 integer, intent(in) :: min_i, min_j, min_k, max_i, max_j, max_k, &
2671 src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
2672 start_x, end_x, start_y, end_y, start_z, end_z, istagger
2673 integer, dimension(src_min_x:src_max_x,src_min_y:src_max_y,2), intent(inout) :: where_maps_to
2674 real, intent(in) :: sr_x, sr_y, maskval, msgval
2675 real, dimension(src_min_x:src_max_x,src_min_y:src_max_y,src_min_z:src_max_z), intent(in) :: tile_array
2676 real, dimension(src_min_x:src_max_x,src_min_y:src_max_y), intent(in), optional :: mask_array
2677 real, dimension(start_x:end_x,start_y:end_y,start_z:end_z), intent(inout) :: dst_array, n
2678 type (bitarray), intent(inout) :: new_pts
2679 character(len=1), intent(in), optional :: mask_relational
2681 ! Local variables
2682 integer :: orig_selected_domain, x_dest, y_dest, i, j, k, center_i, center_j
2683 real :: lat_corner, lon_corner, rx, ry
2685 ! Compute the model grid point that the corners of the rectangle to be
2686 ! processed map to
2687 ! Lower-left corner
2688 if (where_maps_to(min_i,min_j,1) == NOT_PROCESSED) then
2689 orig_selected_domain = iget_selected_domain()
2690 call select_domain(SOURCE_PROJ)
2691 call xytoll(real(min_i), real(min_j), lat_corner, lon_corner, istagger)
2692 call select_domain(1)
2693 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2694 rx = (rx - 1.0)*sr_x + 1.0
2695 ry = (ry - 1.0)*sr_y + 1.0
2696 call select_domain(orig_selected_domain)
2697 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2698 real(start_y) <= ry .and. ry <= real(end_y)) then
2699 where_maps_to(min_i,min_j,1) = nint(rx)
2700 where_maps_to(min_i,min_j,2) = nint(ry)
2701 else
2702 where_maps_to(min_i,min_j,1) = OUTSIDE_DOMAIN
2703 end if
2704 end if
2706 ! Upper-left corner
2707 if (where_maps_to(min_i,max_j,1) == NOT_PROCESSED) then
2708 orig_selected_domain = iget_selected_domain()
2709 call select_domain(SOURCE_PROJ)
2710 call xytoll(real(min_i), real(max_j), lat_corner, lon_corner, istagger)
2711 call select_domain(1)
2712 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2713 rx = (rx - 1.0)*sr_x + 1.0
2714 ry = (ry - 1.0)*sr_y + 1.0
2715 call select_domain(orig_selected_domain)
2716 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2717 real(start_y) <= ry .and. ry <= real(end_y)) then
2718 where_maps_to(min_i,max_j,1) = nint(rx)
2719 where_maps_to(min_i,max_j,2) = nint(ry)
2720 else
2721 where_maps_to(min_i,max_j,1) = OUTSIDE_DOMAIN
2722 end if
2723 end if
2725 ! Upper-right corner
2726 if (where_maps_to(max_i,max_j,1) == NOT_PROCESSED) then
2727 orig_selected_domain = iget_selected_domain()
2728 call select_domain(SOURCE_PROJ)
2729 call xytoll(real(max_i), real(max_j), lat_corner, lon_corner, istagger)
2730 call select_domain(1)
2731 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2732 rx = (rx - 1.0)*sr_x + 1.0
2733 ry = (ry - 1.0)*sr_y + 1.0
2734 call select_domain(orig_selected_domain)
2735 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2736 real(start_y) <= ry .and. ry <= real(end_y)) then
2737 where_maps_to(max_i,max_j,1) = nint(rx)
2738 where_maps_to(max_i,max_j,2) = nint(ry)
2739 else
2740 where_maps_to(max_i,max_j,1) = OUTSIDE_DOMAIN
2741 end if
2742 end if
2744 ! Lower-right corner
2745 if (where_maps_to(max_i,min_j,1) == NOT_PROCESSED) then
2746 orig_selected_domain = iget_selected_domain()
2747 call select_domain(SOURCE_PROJ)
2748 call xytoll(real(max_i), real(min_j), lat_corner, lon_corner, istagger)
2749 call select_domain(1)
2750 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2751 rx = (rx - 1.0)*sr_x + 1.0
2752 ry = (ry - 1.0)*sr_y + 1.0
2753 call select_domain(orig_selected_domain)
2754 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2755 real(start_y) <= ry .and. ry <= real(end_y)) then
2756 where_maps_to(max_i,min_j,1) = nint(rx)
2757 where_maps_to(max_i,min_j,2) = nint(ry)
2758 else
2759 where_maps_to(max_i,min_j,1) = OUTSIDE_DOMAIN
2760 end if
2761 end if
2763 ! If all four corners map to same model grid point, accumulate the
2764 ! entire rectangle
2765 if (where_maps_to(min_i,min_j,1) == where_maps_to(min_i,max_j,1) .and. &
2766 where_maps_to(min_i,min_j,1) == where_maps_to(max_i,max_j,1) .and. &
2767 where_maps_to(min_i,min_j,1) == where_maps_to(max_i,min_j,1) .and. &
2768 where_maps_to(min_i,min_j,2) == where_maps_to(min_i,max_j,2) .and. &
2769 where_maps_to(min_i,min_j,2) == where_maps_to(max_i,max_j,2) .and. &
2770 where_maps_to(min_i,min_j,2) == where_maps_to(max_i,min_j,2) .and. &
2771 where_maps_to(min_i,min_j,1) /= OUTSIDE_DOMAIN) then
2772 x_dest = where_maps_to(min_i,min_j,1)
2773 y_dest = where_maps_to(min_i,min_j,2)
2775 ! If this grid point was already given a value from higher-priority source data,
2776 ! there is nothing to do.
2777 ! if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2779 ! If this grid point has never been given a value by this level of source data,
2780 ! initialize the point
2781 if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2782 do k=min_k,max_k
2783 dst_array(x_dest,y_dest,k) = 0.
2784 end do
2785 end if
2787 ! Sum all the points whose nearest neighbor is this grid point
2788 if (present(mask_array) .and. present(mask_relational)) then
2789 do i=min_i,max_i
2790 do j=min_j,max_j
2791 do k=min_k,max_k
2792 ! Ignore masked/missing values in the source data
2793 if (tile_array(i,j,k) /= msgval) then
2794 if (mask_relational == ' ' .and. mask_array(i,j) /= maskval) then
2795 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2796 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2797 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2798 else if (mask_relational == '<' .and. mask_array(i,j) >= maskval) then
2799 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2800 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2801 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2802 else if (mask_relational == '>' .and. mask_array(i,j) <= maskval) then
2803 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2804 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2805 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2806 end if
2807 end if
2808 end do
2809 end do
2810 end do
2811 else if (present(mask_array)) then
2812 do i=min_i,max_i
2813 do j=min_j,max_j
2814 do k=min_k,max_k
2815 ! Ignore masked/missing values in the source data
2816 if ((tile_array(i,j,k) /= msgval) .and. &
2817 (mask_array(i,j) /= maskval)) then
2818 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2819 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2820 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2821 end if
2822 end do
2823 end do
2824 end do
2825 else
2826 do i=min_i,max_i
2827 do j=min_j,max_j
2828 do k=min_k,max_k
2829 ! Ignore masked/missing values in the source data
2830 if ((tile_array(i,j,k) /= msgval)) then
2831 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2832 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2833 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2834 end if
2835 end do
2836 end do
2837 end do
2838 end if
2840 ! end if
2842 ! Rectangle is a square of four points, and we can simply deal with each of the points
2843 else if (((max_i - min_i + 1) <= 2) .and. ((max_j - min_j + 1) <= 2)) then
2844 do i=min_i,max_i
2845 do j=min_j,max_j
2846 x_dest = where_maps_to(i,j,1)
2847 y_dest = where_maps_to(i,j,2)
2849 if (x_dest /= OUTSIDE_DOMAIN) then
2851 ! if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2852 if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2853 do k=min_k,max_k
2854 dst_array(x_dest,y_dest,k) = 0.
2855 end do
2856 end if
2858 if (present(mask_array) .and. present(mask_relational)) then
2859 do k=min_k,max_k
2860 ! Ignore masked/missing values
2861 if (tile_array(i,j,k) /= msgval) then
2862 if (mask_relational == ' ' .and. mask_array(i,j) /= maskval) then
2863 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2864 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2865 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2866 else if (mask_relational == '<' .and. mask_array(i,j) >= maskval) then
2867 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2868 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2869 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2870 else if (mask_relational == '>' .and. mask_array(i,j) <= maskval) then
2871 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2872 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2873 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2874 end if
2875 end if
2876 end do
2877 else if (present(mask_array)) then
2878 do k=min_k,max_k
2879 ! Ignore masked/missing values
2880 if ((tile_array(i,j,k) /= msgval) .and. &
2881 (mask_array(i,j) /= maskval)) then
2882 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2883 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2884 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2885 end if
2886 end do
2887 else
2888 do k=min_k,max_k
2889 ! Ignore masked/missing values
2890 if ((tile_array(i,j,k) /= msgval)) then
2891 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2892 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2893 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2894 end if
2895 end do
2896 end if
2897 ! end if
2899 end if
2900 end do
2901 end do
2903 ! Not all corners map to the same grid point, and the rectangle contains more than
2904 ! four points
2905 else
2906 center_i = (max_i + min_i)/2
2907 center_j = (max_j + min_j)/2
2909 ! Recursively process lower-left rectangle
2910 call process_continuous_block(tile_array, where_maps_to, &
2911 src_min_x, src_min_y, src_min_z, &
2912 src_max_x, src_max_y, src_max_z, &
2913 min_i, min_j, min_k, &
2914 center_i, center_j, max_k, &
2915 dst_array, n, &
2916 start_x, end_x, start_y, end_y, start_z, end_z, &
2917 istagger, &
2918 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2920 if (center_i < max_i) then
2921 ! Recursively process lower-right rectangle
2922 call process_continuous_block(tile_array, where_maps_to, &
2923 src_min_x, src_min_y, src_min_z, &
2924 src_max_x, src_max_y, src_max_z, &
2925 center_i+1, min_j, min_k, max_i, &
2926 center_j, max_k, &
2927 dst_array, n, &
2928 start_x, end_x, start_y, &
2929 end_y, start_z, end_z, &
2930 istagger, &
2931 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2932 end if
2934 if (center_j < max_j) then
2935 ! Recursively process upper-left rectangle
2936 call process_continuous_block(tile_array, where_maps_to, &
2937 src_min_x, src_min_y, src_min_z, &
2938 src_max_x, src_max_y, src_max_z, &
2939 min_i, center_j+1, min_k, center_i, &
2940 max_j, max_k, &
2941 dst_array, n, &
2942 start_x, end_x, start_y, &
2943 end_y, start_z, end_z, &
2944 istagger, &
2945 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2946 end if
2948 if (center_i < max_i .and. center_j < max_j) then
2949 ! Recursively process upper-right rectangle
2950 call process_continuous_block(tile_array, where_maps_to, &
2951 src_min_x, src_min_y, src_min_z, &
2952 src_max_x, src_max_y, src_max_z, &
2953 center_i+1, center_j+1, min_k, max_i, &
2954 max_j, max_k, &
2955 dst_array, n, &
2956 start_x, end_x, start_y, &
2957 end_y, start_z, end_z, &
2958 istagger, &
2959 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2960 end if
2961 end if
2963 end subroutine process_continuous_block
2965 end module process_domain_module