| blob | 271b8307d68325c4f3469b77a5566d48350d839e |
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
47 ! CWH Initialize local pointer variables
48 nullify(landmask)
49 nullify(xlat)
50 nullify(xlon)
51 nullify(xlat_u)
52 nullify(xlon_u)
53 nullify(xlat_v)
54 nullify(xlon_v)
56 ! Compute number of times that we will process
57 call geth_idts(end_date(n), start_date(n), idiff)
58 call mprintf((idiff < 0),ERROR,'Ending date is earlier than starting date in namelist for domain %i.', i1=n)
60 n_times = idiff / interval_seconds
62 ! Check that the interval evenly divides the range of times to process
63 call mprintf((mod(idiff, interval_seconds) /= 0),WARN, &
64 'In namelist, interval_seconds does not evenly divide '// &
65 '(end_date - start_date) for domain %i. Only %i time periods '// &
66 'will be processed.', i1=n, i2=n_times)
68 ! Initialize the storage module
69 call mprintf(.true.,LOGFILE,'Initializing storage module')
70 call storage_init()
72 !
73 ! Do time-independent processing
74 !
75 call get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, map_proj, &
76 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
77 we_patch_s, we_patch_e, &
78 we_patch_stag_s, we_patch_stag_e, &
79 sn_patch_s, sn_patch_e, &
80 sn_patch_stag_s, sn_patch_stag_e, &
81 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
82 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
83 mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
84 grid_id, parent_id, &
85 i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
86 parent_grid_ratio, sub_x, sub_y, &
87 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
88 pole_lat, pole_lon, dom_dx, dom_dy, landmask, xlat, xlon, xlat_u, xlon_u, &
89 xlat_v, xlon_v, corner_lats, corner_lons, title)
92 allocate(output_flags(num_entries))
93 allocate(got_const_field(num_entries))
95 do i=1,num_entries
96 output_flags(i) = ' '
97 got_const_field(i) = .false.
98 end do
100 ! This call is to process the constant met fields (SST or SEAICE, for example)
101 ! That we process constant fields is indicated by the first argument
102 call process_single_met_time(.true., temp_date, n, extra_row, extra_col, xlat, xlon, &
103 xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
104 title, dyn_opt, &
105 west_east_dim, south_north_dim, &
106 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
107 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
108 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
109 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
110 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
111 got_const_field, &
112 map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
113 grid_id, parent_id, i_parent_start, &
114 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
115 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
116 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
117 corner_lats, corner_lons, output_flags, 0)
119 !
120 ! Begin time-dependent processing
121 !
123 allocate(td_output_flags(num_entries))
124 allocate(got_this_field (num_entries))
126 ! Loop over all times to be processed for this domain
127 do t=0,n_times
129 call geth_newdate(valid_date, trim(start_date(n)), t*interval_seconds)
130 temp_date = ' '
132 if (mod(interval_seconds,3600) == 0) then
133 write(temp_date,'(a13)') valid_date(1:10)//'_'//valid_date(12:13)
134 else if (mod(interval_seconds,60) == 0) then
135 write(temp_date,'(a16)') valid_date(1:10)//'_'//valid_date(12:16)
136 else
137 write(temp_date,'(a19)') valid_date(1:10)//'_'//valid_date(12:19)
138 end if
140 call mprintf(.true.,STDOUT, ' Processing %s', s1=trim(temp_date))
141 call mprintf(.true.,LOGFILE, 'Preparing to process output time %s', s1=temp_date)
143 do i=1,num_entries
144 td_output_flags(i) = output_flags(i)
145 got_this_field(i) = got_const_field(i)
146 end do
148 if (t > 0) then
149 process_bdy_width = process_only_bdy
150 else
151 process_bdy_width = 0
152 end if
154 call process_single_met_time(.false., temp_date, n, extra_row, extra_col, xlat, xlon, &
155 xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
156 title, dyn_opt, &
157 west_east_dim, south_north_dim, &
158 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
159 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
160 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
161 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
162 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
163 got_this_field, &
164 map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
165 grid_id, parent_id, i_parent_start, &
166 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
167 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
168 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
169 corner_lats, corner_lons, td_output_flags, process_bdy_width)
171 end do ! Loop over n_times
174 deallocate(td_output_flags)
175 deallocate(got_this_field)
177 deallocate(output_flags)
178 deallocate(got_const_field)
180 call storage_delete_all()
182 end subroutine process_domain
185 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
186 ! Name: get_static_fields
187 !
188 ! Purpose:
189 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
190 subroutine get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, &
191 map_proj, &
192 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
193 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
194 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
195 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
196 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
197 mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
198 grid_id, parent_id, &
199 i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
200 parent_grid_ratio, sub_x, sub_y, &
201 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
202 pole_lat, pole_lon, dom_dx, dom_dy, landmask, xlat, xlon, xlat_u, xlon_u, &
203 xlat_v, xlon_v, corner_lats, corner_lons, title)
205 use gridinfo_module
206 use input_module
207 use llxy_module
208 use parallel_module
209 use storage_module
210 use module_debug
212 implicit none
214 ! Arguments
215 integer, intent(in) :: n
216 integer, intent(inout) :: dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, &
217 map_proj, &
218 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
219 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
220 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
221 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
222 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
223 is_water, is_lake, is_ice, is_urban, i_soilwater, grid_id, parent_id, &
224 i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
225 parent_grid_ratio, sub_x, sub_y, num_land_cat
226 real, pointer, dimension(:,:) :: landmask
227 real, intent(inout) :: cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
228 dom_dx, dom_dy, pole_lat, pole_lon
229 real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
230 real, dimension(16), intent(out) :: corner_lats, corner_lons
231 character (len=128), intent(inout) :: title, mminlu
233 ! Local variables
234 integer :: istatus, i, j, k, sp1, ep1, sp2, ep2, sp3, ep3, &
235 lh_mult, rh_mult, bh_mult, th_mult, subx, suby
236 integer :: we_mem_subgrid_s, we_mem_subgrid_e, &
237 sn_mem_subgrid_s, sn_mem_subgrid_e
238 integer :: we_patch_subgrid_s, we_patch_subgrid_e, &
239 sn_patch_subgrid_s, sn_patch_subgrid_e
240 real, pointer, dimension(:,:,:) :: real_array
241 character (len=3) :: memorder
242 character (len=128) :: grid_type, datestr, cname, stagger, cunits, cdesc
243 character (len=128), dimension(3) :: dimnames
244 type (fg_input) :: field
246 ! CWH Initialize local pointer variables
247 nullify(real_array)
249 ! Initialize the input module to read static input data for this domain
250 call mprintf(.true.,LOGFILE,'Opening static input file.')
251 call input_init(n, istatus)
252 call mprintf((istatus /= 0),ERROR, 'input_init(): Error opening input for domain %i.', i1=n)
254 ! Read global attributes from the static data input file
255 call mprintf(.true.,LOGFILE,'Reading static global attributes.')
256 call read_global_attrs(title, datestr, grid_type, dyn_opt, west_east_dim, &
257 south_north_dim, bottom_top_dim, &
258 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
259 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
260 map_proj, mminlu, num_land_cat, &
261 is_water, is_lake, is_ice, is_urban, i_soilwater, &
262 grid_id, parent_id, i_parent_start, &
263 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
264 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, &
265 truelat2, pole_lat, pole_lon, parent_grid_ratio, &
266 corner_lats, corner_lons, sub_x, sub_y)
268 we_dom_s = 1
269 sn_dom_s = 1
270 if (grid_type(1:1) == 'C') then
271 we_dom_e = west_east_dim - 1
272 sn_dom_e = south_north_dim - 1
273 else if (grid_type(1:1) == 'E') then
274 we_dom_e = west_east_dim
275 sn_dom_e = south_north_dim
276 end if
278 ! Given the full dimensions of this domain, find out the range of indices
279 ! that will be worked on by this processor. This information is given by
280 ! my_minx, my_miny, my_maxx, my_maxy
281 call parallel_get_tile_dims(west_east_dim, south_north_dim)
283 ! Must figure out patch dimensions from info in parallel module
284 if (nprocs > 1 .and. .not. do_tiled_input) then
286 we_patch_s = my_minx
287 we_patch_stag_s = my_minx
288 we_patch_e = my_maxx - 1
289 sn_patch_s = my_miny
290 sn_patch_stag_s = my_miny
291 sn_patch_e = my_maxy - 1
293 if (gridtype == 'C') then
294 if (my_x /= nproc_x - 1) then
295 we_patch_e = we_patch_e + 1
296 we_patch_stag_e = we_patch_e
297 else
298 we_patch_stag_e = we_patch_e + 1
299 end if
300 if (my_y /= nproc_y - 1) then
301 sn_patch_e = sn_patch_e + 1
302 sn_patch_stag_e = sn_patch_e
303 else
304 sn_patch_stag_e = sn_patch_e + 1
305 end if
306 else if (gridtype == 'E') then
307 we_patch_e = we_patch_e + 1
308 sn_patch_e = sn_patch_e + 1
309 we_patch_stag_e = we_patch_e
310 sn_patch_stag_e = sn_patch_e
311 end if
313 end if
315 ! Compute multipliers for halo width; these must be 0/1
316 if (my_x /= 0) then
317 lh_mult = 1
318 else
319 lh_mult = 0
320 end if
321 if (my_x /= (nproc_x-1)) then
322 rh_mult = 1
323 else
324 rh_mult = 0
325 end if
326 if (my_y /= 0) then
327 bh_mult = 1
328 else
329 bh_mult = 0
330 end if
331 if (my_y /= (nproc_y-1)) then
332 th_mult = 1
333 else
334 th_mult = 0
335 end if
337 we_mem_s = we_patch_s - HALO_WIDTH*lh_mult
338 we_mem_e = we_patch_e + HALO_WIDTH*rh_mult
339 sn_mem_s = sn_patch_s - HALO_WIDTH*bh_mult
340 sn_mem_e = sn_patch_e + HALO_WIDTH*th_mult
341 we_mem_stag_s = we_patch_stag_s - HALO_WIDTH*lh_mult
342 we_mem_stag_e = we_patch_stag_e + HALO_WIDTH*rh_mult
343 sn_mem_stag_s = sn_patch_stag_s - HALO_WIDTH*bh_mult
344 sn_mem_stag_e = sn_patch_stag_e + HALO_WIDTH*th_mult
346 ! Initialize a proj_info type for the destination grid projection. This will
347 ! primarily be used for rotating Earth-relative winds to grid-relative winds
348 call set_domain_projection(map_proj, stand_lon, truelat1, truelat2, &
349 dom_dx, dom_dy, dom_dx, dom_dy, west_east_dim, &
350 south_north_dim, real(west_east_dim)/2., &
351 real(south_north_dim)/2.,cen_lat, cen_lon, &
352 cen_lat, cen_lon)
354 ! Read static fields using the input module; we know that there are no more
355 ! fields to be read when read_next_field() returns a non-zero status.
356 istatus = 0
357 do while (istatus == 0)
358 call read_next_field(sp1, ep1, sp2, ep2, sp3, ep3, cname, cunits, cdesc, &
359 memorder, stagger, dimnames, subx, suby, &
360 real_array, istatus)
361 if (istatus == 0) then
363 call mprintf(.true.,LOGFILE, 'Read in static field %s.',s1=cname)
365 ! We will also keep copies in core of the lat/lon arrays, for use in
366 ! interpolation of the met fields to the model grid.
367 ! For now, we assume that the lat/lon arrays will have known field names
368 if (index(cname, 'XLAT_M') /= 0 .and. &
369 len_trim(cname) == len_trim('XLAT_M')) then
370 allocate(xlat(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
371 xlat(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
372 call exchange_halo_r(xlat, &
373 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
374 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
376 else if (index(cname, 'XLONG_M') /= 0 .and. &
377 len_trim(cname) == len_trim('XLONG_M')) then
378 allocate(xlon(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
379 xlon(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
380 call exchange_halo_r(xlon, &
381 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
382 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
384 else if (index(cname, 'XLAT_U') /= 0 .and. &
385 len_trim(cname) == len_trim('XLAT_U')) then
386 allocate(xlat_u(we_mem_stag_s:we_mem_stag_e,sn_mem_s:sn_mem_e))
387 xlat_u(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
388 call exchange_halo_r(xlat_u, &
389 we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
390 we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
392 else if (index(cname, 'XLONG_U') /= 0 .and. &
393 len_trim(cname) == len_trim('XLONG_U')) then
394 allocate(xlon_u(we_mem_stag_s:we_mem_stag_e,sn_mem_s:sn_mem_e))
395 xlon_u(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
396 call exchange_halo_r(xlon_u, &
397 we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
398 we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
400 else if (index(cname, 'XLAT_V') /= 0 .and. &
401 len_trim(cname) == len_trim('XLAT_V')) then
402 allocate(xlat_v(we_mem_s:we_mem_e,sn_mem_stag_s:sn_mem_stag_e))
403 xlat_v(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(:,:,1)
404 call exchange_halo_r(xlat_v, &
405 we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
406 we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
408 else if (index(cname, 'XLONG_V') /= 0 .and. &
409 len_trim(cname) == len_trim('XLONG_V')) then
410 allocate(xlon_v(we_mem_s:we_mem_e,sn_mem_stag_s:sn_mem_stag_e))
411 xlon_v(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(:,:,1)
412 call exchange_halo_r(xlon_v, &
413 we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
414 we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
416 else if (index(cname, 'LANDMASK') /= 0 .and. &
417 len_trim(cname) == len_trim('LANDMASK')) then
418 allocate(landmask(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
419 landmask(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
420 call exchange_halo_r(landmask, &
421 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
422 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
424 end if
426 if (subx > 1) then
427 we_mem_subgrid_s = (we_mem_s + HALO_WIDTH*lh_mult - 1)*subx - HALO_WIDTH*lh_mult + 1
428 we_mem_subgrid_e = (we_mem_e + (1-rh_mult) - HALO_WIDTH*rh_mult )*subx + HALO_WIDTH*rh_mult
429 we_patch_subgrid_s = (we_patch_s - 1)*subx + 1
430 we_patch_subgrid_e = (we_patch_e + (1-rh_mult) )*subx
431 end if
432 if (suby > 1) then
433 sn_mem_subgrid_s = (sn_mem_s + HALO_WIDTH*bh_mult - 1)*suby - HALO_WIDTH*bh_mult + 1
434 sn_mem_subgrid_e = (sn_mem_e + (1-th_mult) - HALO_WIDTH*th_mult )*suby + HALO_WIDTH*th_mult
435 sn_patch_subgrid_s = (sn_patch_s - 1)*suby + 1
436 sn_patch_subgrid_e = (sn_patch_e + (1-th_mult) )*suby
437 end if
439 ! Having read in a field, we write each level individually to the
440 ! storage module; levels will be reassembled later on when they
441 ! are written.
442 do k=sp3,ep3
443 field%header%sr_x=subx
444 field%header%sr_y=suby
445 field%header%version = 1
446 field%header%date = start_date(n)
447 field%header%time_dependent = .false.
448 field%header%mask_field = .false.
449 field%header%forecast_hour = 0.0
450 field%header%fg_source = 'geogrid_model'
451 field%header%field = cname
452 field%header%units = cunits
453 field%header%description = cdesc
454 field%header%vertical_coord = dimnames(3)
455 field%header%vertical_level = k
456 field%header%array_order = memorder
457 field%header%is_wind_grid_rel = .true.
458 field%header%array_has_missing_values = .false.
459 if (gridtype == 'C') then
460 if (subx > 1 .or. suby > 1) then
461 field%map%stagger = M
462 field%header%dim1(1) = we_mem_subgrid_s
463 field%header%dim1(2) = we_mem_subgrid_e
464 field%header%dim2(1) = sn_mem_subgrid_s
465 field%header%dim2(2) = sn_mem_subgrid_e
466 else if (trim(stagger) == 'M') then
467 field%map%stagger = M
468 field%header%dim1(1) = we_mem_s
469 field%header%dim1(2) = we_mem_e
470 field%header%dim2(1) = sn_mem_s
471 field%header%dim2(2) = sn_mem_e
472 else if (trim(stagger) == 'U') then
473 field%map%stagger = U
474 field%header%dim1(1) = we_mem_stag_s
475 field%header%dim1(2) = we_mem_stag_e
476 field%header%dim2(1) = sn_mem_s
477 field%header%dim2(2) = sn_mem_e
478 else if (trim(stagger) == 'V') then
479 field%map%stagger = V
480 field%header%dim1(1) = we_mem_s
481 field%header%dim1(2) = we_mem_e
482 field%header%dim2(1) = sn_mem_stag_s
483 field%header%dim2(2) = sn_mem_stag_e
484 end if
485 else if (gridtype == 'E') then
486 if (trim(stagger) == 'M') then
487 field%map%stagger = HH
488 else if (trim(stagger) == 'V') then
489 field%map%stagger = VV
490 end if
491 field%header%dim1(1) = we_mem_s
492 field%header%dim1(2) = we_mem_e
493 field%header%dim2(1) = sn_mem_s
494 field%header%dim2(2) = sn_mem_e
495 end if
497 allocate(field%valid_mask)
499 if (subx > 1 .or. suby > 1) then
500 allocate(field%r_arr(we_mem_subgrid_s:we_mem_subgrid_e,&
501 sn_mem_subgrid_s:sn_mem_subgrid_e))
502 field%r_arr(we_patch_subgrid_s:we_patch_subgrid_e,sn_patch_subgrid_s:sn_patch_subgrid_e) = &
503 real_array(sp1:ep1,sp2:ep2,k)
504 call exchange_halo_r(field%r_arr, &
505 we_mem_subgrid_s, we_mem_subgrid_e, sn_mem_subgrid_s, sn_mem_subgrid_e, 1, 1, &
506 we_patch_subgrid_s, we_patch_subgrid_e, sn_patch_subgrid_s, sn_patch_subgrid_e, 1, 1)
507 call bitarray_create(field%valid_mask, &
508 (we_mem_subgrid_e-we_mem_subgrid_s)+1, &
509 (sn_mem_subgrid_e-sn_mem_subgrid_s)+1)
510 do j=1,(sn_mem_subgrid_e-sn_mem_subgrid_s)+1
511 do i=1,(we_mem_subgrid_e-we_mem_subgrid_s)+1
512 call bitarray_set(field%valid_mask, i, j)
513 end do
514 end do
516 else if (field%map%stagger == M .or. &
517 field%map%stagger == HH .or. &
518 field%map%stagger == VV) then
519 allocate(field%r_arr(we_mem_s:we_mem_e,&
520 sn_mem_s:sn_mem_e))
521 field%r_arr(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(sp1:ep1,sp2:ep2,k)
522 call exchange_halo_r(field%r_arr, &
523 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
524 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
525 call bitarray_create(field%valid_mask, &
526 (we_mem_e-we_mem_s)+1, &
527 (sn_mem_e-sn_mem_s)+1)
528 do j=1,(sn_mem_e-sn_mem_s)+1
529 do i=1,(we_mem_e-we_mem_s)+1
530 call bitarray_set(field%valid_mask, i, j)
531 end do
532 end do
533 else if (field%map%stagger == U) then
534 allocate(field%r_arr(we_mem_stag_s:we_mem_stag_e,&
535 sn_mem_s:sn_mem_e))
536 field%r_arr(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(sp1:ep1,sp2:ep2,k)
537 call exchange_halo_r(field%r_arr, &
538 we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
539 we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
540 call bitarray_create(field%valid_mask, &
541 (we_mem_stag_e-we_mem_stag_s)+1, &
542 (sn_mem_e-sn_mem_s)+1)
543 do j=1,(sn_mem_e-sn_mem_s)+1
544 do i=1,(we_mem_stag_e-we_mem_stag_s)+1
545 call bitarray_set(field%valid_mask, i, j)
546 end do
547 end do
548 else if (field%map%stagger == V) then
549 allocate(field%r_arr(we_mem_s:we_mem_e,&
550 sn_mem_stag_s:sn_mem_stag_e))
551 field%r_arr(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(sp1:ep1,sp2:ep2,k)
552 call exchange_halo_r(field%r_arr, &
553 we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
554 we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
555 call bitarray_create(field%valid_mask, &
556 (we_mem_e-we_mem_s)+1, &
557 (sn_mem_stag_e-sn_mem_stag_s)+1)
558 do j=1,(sn_mem_stag_e-sn_mem_stag_s)+1
559 do i=1,(we_mem_e-we_mem_s)+1
560 call bitarray_set(field%valid_mask, i, j)
561 end do
562 end do
563 end if
565 nullify(field%modified_mask)
567 call storage_put_field(field)
569 end do
571 end if
572 end do
574 ! Done reading all static fields for this domain
575 call input_close()
577 end subroutine get_static_fields
580 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
581 ! Name: process_single_met_time
582 !
583 ! Purpose:
584 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
585 subroutine process_single_met_time(do_const_processing, &
586 temp_date, n, extra_row, extra_col, xlat, xlon, &
587 xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
588 title, dyn_opt, &
589 west_east_dim, south_north_dim, &
590 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
591 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
592 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
593 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
594 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
595 got_this_field, &
596 map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
597 grid_id, parent_id, i_parent_start, &
598 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
599 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
600 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
601 corner_lats, corner_lons, output_flags, process_bdy_width)
603 use bitarray_module
604 use gridinfo_module
605 use interp_module
606 use interp_option_module
607 use llxy_module
608 use misc_definitions_module
609 use module_debug
610 use output_module
611 use parallel_module
612 use read_met_module
613 use rotate_winds_module
614 use storage_module
616 implicit none
618 ! Arguments
619 logical, intent(in) :: do_const_processing
620 integer, intent(in) :: n, dyn_opt, west_east_dim, south_north_dim, map_proj, &
621 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
622 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
623 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
624 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
625 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
626 is_water, is_lake, is_ice, is_urban, i_soilwater, &
627 grid_id, parent_id, i_parent_start, j_parent_start, &
628 i_parent_end, j_parent_end, parent_grid_ratio, sub_x, sub_y, num_land_cat, &
629 process_bdy_width
630 ! BUG: Should we be passing these around as pointers, or just declare them as arrays?
631 real, pointer, dimension(:,:) :: landmask
632 real, intent(in) :: dom_dx, dom_dy, cen_lat, moad_cen_lat, cen_lon, stand_lon, &
633 truelat1, truelat2, pole_lat, pole_lon
634 real, dimension(16), intent(in) :: corner_lats, corner_lons
635 real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
636 logical, intent(in) :: extra_row, extra_col
637 logical, dimension(:), intent(inout) :: got_this_field
638 character (len=19), intent(in) :: temp_date
639 character (len=128), intent(in) :: mminlu
640 character (len=128), dimension(:), intent(inout) :: output_flags
642 ! BUG: Move this constant to misc_definitions_module?
643 integer, parameter :: BDR_WIDTH = 3
645 ! Local variables
646 integer :: istatus, iqstatus, fg_idx, idx, idxt, i, j, bottom_top_dim, &
647 sm1, em1, sm2, em2, sm3, em3, &
648 sp1, ep1, sp2, ep2, sp3, ep3, &
649 sd1, ed1, sd2, ed2, sd3, ed3, &
650 u_idx, bdr_wdth
651 integer :: nmet_flags
652 integer :: num_metgrid_soil_levs
653 integer, pointer, dimension(:) :: soil_levels
654 real :: rx, ry
655 real :: threshold
656 logical :: do_gcell_interp
657 integer, pointer, dimension(:) :: u_levels, v_levels
658 real, pointer, dimension(:,:) :: halo_slab
659 real, pointer, dimension(:,:,:) :: real_array
660 character (len=19) :: output_date
661 character (len=128) :: cname, title
662 character (len=MAX_FILENAME_LEN) :: input_name
663 character (len=128), allocatable, dimension(:) :: met_flags
664 type (fg_input) :: field, u_field, v_field
665 type (met_data) :: fg_data
667 ! CWH Initialize local pointer variables
668 nullify(soil_levels)
669 nullify(u_levels)
670 nullify(v_levels)
671 nullify(halo_slab)
672 nullify(real_array)
675 ! For this time, we need to process all first-guess filename roots. When we
676 ! hit a root containing a '*', we assume we have hit the end of the list
677 fg_idx = 1
678 if (do_const_processing) then
679 input_name = constants_name(fg_idx)
680 else
681 input_name = fg_name(fg_idx)
682 end if
683 do while (input_name /= '*')
685 call mprintf(.true.,STDOUT, ' %s', s1=input_name)
686 call mprintf(.true.,LOGFILE, 'Getting input fields from %s', s1=input_name)
688 ! Do a first pass through this fg source to get all mask fields used
689 ! during interpolation
690 call get_interp_masks(trim(input_name), do_const_processing, temp_date)
692 istatus = 0
694 ! Initialize the module for reading in the met fields
695 call read_met_init(trim(input_name), do_const_processing, temp_date, istatus)
697 if (istatus == 0) then
699 ! Process all fields and levels from the current file; read_next_met_field()
700 ! will return a non-zero status when there are no more fields to be read.
701 do while (istatus == 0)
703 call read_next_met_field(fg_data, istatus)
705 if (istatus == 0) then
707 ! Find index into fieldname, interp_method, masked, and fill_missing
708 ! of the current field
709 idxt = num_entries + 1
710 do idx=1,num_entries
711 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
712 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
714 got_this_field(idx) = .true.
716 if (index(input_name,trim(from_input(idx))) /= 0 .or. &
717 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
718 idxt = idx
719 end if
721 end if
722 end do
723 idx = idxt
724 if (idx > num_entries) idx = num_entries ! The last entry is a default
726 ! Do we need to rename this field?
727 if (output_name(idx) /= ' ') then
728 fg_data%field = output_name(idx)(1:9)
730 idxt = num_entries + 1
731 do idx=1,num_entries
732 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
733 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
735 got_this_field(idx) = .true.
737 if (index(input_name,trim(from_input(idx))) /= 0 .or. &
738 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
739 idxt = idx
740 end if
742 end if
743 end do
744 idx = idxt
745 if (idx > num_entries) idx = num_entries ! The last entry is a default
746 end if
748 ! Do a simple check to see whether this is a global dataset
749 ! Note that we do not currently support regional Gaussian grids
750 if ((fg_data%iproj == PROJ_LATLON .and. abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) &
751 .or. (fg_data%iproj == PROJ_GAUSS)) then
752 bdr_wdth = BDR_WIDTH
753 allocate(halo_slab(1-BDR_WIDTH:fg_data%nx+BDR_WIDTH,1:fg_data%ny))
755 halo_slab(1:fg_data%nx, 1:fg_data%ny) = &
756 fg_data%slab(1:fg_data%nx, 1:fg_data%ny)
758 halo_slab(1-BDR_WIDTH:0, 1:fg_data%ny) = &
759 fg_data%slab(fg_data%nx-BDR_WIDTH+1:fg_data%nx, 1:fg_data%ny)
761 halo_slab(fg_data%nx+1:fg_data%nx+BDR_WIDTH, 1:fg_data%ny) = &
762 fg_data%slab(1:BDR_WIDTH, 1:fg_data%ny)
764 deallocate(fg_data%slab)
765 else
766 bdr_wdth = 0
767 halo_slab => fg_data%slab
768 nullify(fg_data%slab)
769 end if
771 call mprintf(.true.,LOGFILE,'Processing %s at level %f.',s1=fg_data%field,f1=fg_data%xlvl)
773 call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, &
774 fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, &
775 fg_data%deltalon, fg_data%starti, fg_data%startj, &
776 fg_data%startlat, fg_data%startlon, earth_radius=fg_data%earth_radius*1000.)
778 ! Initialize fg_input structure to store the field
779 field%header%version = 1
780 field%header%date = fg_data%hdate//' '
781 if (do_const_processing) then
782 field%header%time_dependent = .false.
783 else
784 field%header%time_dependent = .true.
785 end if
786 field%header%forecast_hour = fg_data%xfcst
787 field%header%fg_source = 'FG'
788 field%header%field = ' '
789 field%header%field(1:9) = fg_data%field
790 field%header%units = ' '
791 field%header%units(1:25) = fg_data%units
792 field%header%description = ' '
793 field%header%description(1:46) = fg_data%desc
794 call get_z_dim_name(fg_data%field,field%header%vertical_coord)
795 field%header%vertical_level = nint(fg_data%xlvl)
796 field%header%sr_x = 1
797 field%header%sr_y = 1
798 field%header%array_order = 'XY '
799 field%header%is_wind_grid_rel = fg_data%is_wind_grid_rel
800 field%header%array_has_missing_values = .false.
801 nullify(field%r_arr)
802 nullify(field%valid_mask)
803 nullify(field%modified_mask)
805 if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
806 output_flags(idx) = flag_in_output(idx)
807 end if
809 ! If we should not output this field, just list it as a mask field
810 if (output_this_field(idx)) then
811 field%header%mask_field = .false.
812 else
813 field%header%mask_field = .true.
814 end if
816 !
817 ! Before actually doing any interpolation to the model grid, we must check
818 ! whether we will be using the average_gcell interpolator that averages all
819 ! source points in each model grid cell
820 !
821 do_gcell_interp = .false.
822 if (index(interp_method(idx),'average_gcell') /= 0) then
824 call get_gcell_threshold(interp_method(idx), threshold, istatus)
825 if (istatus == 0) then
826 if (fg_data%dx == 0. .and. fg_data%dy == 0. .and. &
827 fg_data%deltalat /= 0. .and. fg_data%deltalon /= 0.) then
828 fg_data%dx = abs(fg_data%deltalon)
829 fg_data%dy = abs(fg_data%deltalat)
830 else
831 ! BUG: Need to more correctly handle dx/dy in meters.
832 fg_data%dx = fg_data%dx / 111000. ! Convert meters to approximate degrees
833 fg_data%dy = fg_data%dy / 111000.
834 end if
835 if (gridtype == 'C') then
836 if (threshold*max(fg_data%dx,fg_data%dy)*111. <= max(dom_dx,dom_dy)/1000.) &
837 do_gcell_interp = .true.
838 else if (gridtype == 'E') then
839 if (threshold*max(fg_data%dx,fg_data%dy) <= max(dom_dx,dom_dy)) &
840 do_gcell_interp = .true.
841 end if
842 end if
843 end if
845 ! Interpolate to U staggering
846 if (output_stagger(idx) == U) then
848 call storage_query_field(field, iqstatus)
849 if (iqstatus == 0) then
850 call storage_get_field(field, iqstatus)
851 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
852 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
853 if (associated(field%modified_mask)) then
854 call bitarray_destroy(field%modified_mask)
855 nullify(field%modified_mask)
856 end if
857 else
858 allocate(field%valid_mask)
859 call bitarray_create(field%valid_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
860 end if
862 ! Save a copy of the fg_input structure for the U field so that we can find it later
863 if (is_u_field(idx)) call dup(field, u_field)
865 allocate(field%modified_mask)
866 call bitarray_create(field%modified_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
868 if (do_const_processing .or. field%header%time_dependent) then
869 call interp_met_field(input_name, fg_data%field, U, M, &
870 field, xlat_u, xlon_u, we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, &
871 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
872 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
873 field%modified_mask, process_bdy_width)
874 else
875 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
876 end if
878 ! Interpolate to V staggering
879 else if (output_stagger(idx) == V) then
881 call storage_query_field(field, iqstatus)
882 if (iqstatus == 0) then
883 call storage_get_field(field, iqstatus)
884 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
885 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
886 if (associated(field%modified_mask)) then
887 call bitarray_destroy(field%modified_mask)
888 nullify(field%modified_mask)
889 end if
890 else
891 allocate(field%valid_mask)
892 call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
893 end if
895 ! Save a copy of the fg_input structure for the V field so that we can find it later
896 if (is_v_field(idx)) call dup(field, v_field)
898 allocate(field%modified_mask)
899 call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
901 if (do_const_processing .or. field%header%time_dependent) then
902 call interp_met_field(input_name, fg_data%field, V, M, &
903 field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
904 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
905 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
906 field%modified_mask, process_bdy_width)
907 else
908 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
909 end if
911 ! Interpolate to VV staggering
912 else if (output_stagger(idx) == VV) then
914 call storage_query_field(field, iqstatus)
915 if (iqstatus == 0) then
916 call storage_get_field(field, iqstatus)
917 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
918 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
919 if (associated(field%modified_mask)) then
920 call bitarray_destroy(field%modified_mask)
921 nullify(field%modified_mask)
922 end if
923 else
924 allocate(field%valid_mask)
925 call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
926 end if
928 ! Save a copy of the fg_input structure for the U field so that we can find it later
929 if (is_u_field(idx)) call dup(field, u_field)
931 ! Save a copy of the fg_input structure for the V field so that we can find it later
932 if (is_v_field(idx)) call dup(field, v_field)
934 allocate(field%modified_mask)
935 call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
937 if (do_const_processing .or. field%header%time_dependent) then
938 call interp_met_field(input_name, fg_data%field, VV, M, &
939 field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
940 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
941 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
942 field%modified_mask, process_bdy_width)
943 else
944 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
945 end if
947 ! All other fields interpolated to M staggering for C grid, H staggering for E grid
948 else
950 call storage_query_field(field, iqstatus)
951 if (iqstatus == 0) then
952 call storage_get_field(field, iqstatus)
953 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
954 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
955 if (associated(field%modified_mask)) then
956 call bitarray_destroy(field%modified_mask)
957 nullify(field%modified_mask)
958 end if
959 else
960 allocate(field%valid_mask)
961 call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
962 end if
964 allocate(field%modified_mask)
965 call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
967 if (do_const_processing .or. field%header%time_dependent) then
968 if (gridtype == 'C') then
969 call interp_met_field(input_name, fg_data%field, M, M, &
970 field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
971 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
972 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
973 field%modified_mask, process_bdy_width, landmask)
975 else if (gridtype == 'E') then
976 call interp_met_field(input_name, fg_data%field, HH, M, &
977 field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
978 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
979 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
980 field%modified_mask, process_bdy_width, landmask)
981 end if
982 else
983 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
984 end if
986 end if
988 call bitarray_merge(field%valid_mask, field%modified_mask)
990 deallocate(halo_slab)
992 ! Store the interpolated field
993 call storage_put_field(field)
995 call pop_source_projection()
997 end if
998 end do
1000 call read_met_close()
1002 call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, &
1003 fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, &
1004 fg_data%deltalon, fg_data%starti, fg_data%startj, &
1005 fg_data%startlat, fg_data%startlon, earth_radius=fg_data%earth_radius*1000.)
1007 !
1008 ! If necessary, rotate winds to earth-relative for this fg source
1009 !
1011 call storage_get_levels(u_field, u_levels)
1012 call storage_get_levels(v_field, v_levels)
1014 if (associated(u_levels) .and. associated(v_levels)) then
1015 u_idx = 1
1016 do u_idx = 1, size(u_levels)
1017 u_field%header%vertical_level = u_levels(u_idx)
1018 call storage_get_field(u_field, istatus)
1019 v_field%header%vertical_level = v_levels(u_idx)
1020 call storage_get_field(v_field, istatus)
1022 if (associated(u_field%modified_mask) .and. &
1023 associated(v_field%modified_mask)) then
1025 if (u_field%header%is_wind_grid_rel) then
1026 if (gridtype == 'C') then
1027 call map_to_met(u_field%r_arr, u_field%modified_mask, &
1028 v_field%r_arr, v_field%modified_mask, &
1029 we_mem_stag_s, sn_mem_s, &
1030 we_mem_stag_e, sn_mem_e, &
1031 we_mem_s, sn_mem_stag_s, &
1032 we_mem_e, sn_mem_stag_e, &
1033 xlon_u, xlon_v, xlat_u, xlat_v)
1034 else if (gridtype == 'E') then
1035 call map_to_met_nmm(u_field%r_arr, u_field%modified_mask, &
1036 v_field%r_arr, v_field%modified_mask, &
1037 we_mem_s, sn_mem_s, &
1038 we_mem_e, sn_mem_e, &
1039 xlat_v, xlon_v)
1040 end if
1041 end if
1043 call bitarray_destroy(u_field%modified_mask)
1044 call bitarray_destroy(v_field%modified_mask)
1045 nullify(u_field%modified_mask)
1046 nullify(v_field%modified_mask)
1047 call storage_put_field(u_field)
1048 call storage_put_field(v_field)
1049 end if
1051 end do
1053 deallocate(u_levels)
1054 deallocate(v_levels)
1056 end if
1058 call pop_source_projection()
1060 else
1061 if (do_const_processing) then
1062 call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=input_name)
1063 else
1064 call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=trim(input_name)//':'//trim(temp_date))
1065 end if
1066 end if
1068 fg_idx = fg_idx + 1
1069 if (do_const_processing) then
1070 input_name = constants_name(fg_idx)
1071 else
1072 input_name = fg_name(fg_idx)
1073 end if
1074 end do ! while (input_name /= '*')
1076 !
1077 ! Rotate winds from earth-relative to grid-relative
1078 !
1080 call storage_get_levels(u_field, u_levels)
1081 call storage_get_levels(v_field, v_levels)
1083 if (associated(u_levels) .and. associated(v_levels)) then
1084 u_idx = 1
1085 do u_idx = 1, size(u_levels)
1086 u_field%header%vertical_level = u_levels(u_idx)
1087 call storage_get_field(u_field, istatus)
1088 v_field%header%vertical_level = v_levels(u_idx)
1089 call storage_get_field(v_field, istatus)
1091 if (gridtype == 'C') then
1092 call met_to_map(u_field%r_arr, u_field%valid_mask, &
1093 v_field%r_arr, v_field%valid_mask, &
1094 we_mem_stag_s, sn_mem_s, &
1095 we_mem_stag_e, sn_mem_e, &
1096 we_mem_s, sn_mem_stag_s, &
1097 we_mem_e, sn_mem_stag_e, &
1098 xlon_u, xlon_v, xlat_u, xlat_v)
1099 else if (gridtype == 'E') then
1100 call met_to_map_nmm(u_field%r_arr, u_field%valid_mask, &
1101 v_field%r_arr, v_field%valid_mask, &
1102 we_mem_s, sn_mem_s, &
1103 we_mem_e, sn_mem_e, &
1104 xlat_v, xlon_v)
1105 end if
1107 end do
1109 deallocate(u_levels)
1110 deallocate(v_levels)
1112 end if
1114 if (do_const_processing) return
1116 !
1117 ! Now that we have all degribbed fields, we build a 3-d pressure field, and fill in any
1118 ! missing levels in the other 3-d fields
1119 !
1120 call mprintf(.true.,LOGFILE,'Filling missing levels.')
1121 call fill_missing_levels(output_flags)
1123 call mprintf(.true.,LOGFILE,'Creating derived fields.')
1124 call create_derived_fields(gridtype, fg_data%hdate, fg_data%xfcst, &
1125 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1126 we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e, &
1127 got_this_field, output_flags)
1129 !
1130 ! Check that every mandatory field was found in input data
1131 !
1132 do i=1,num_entries
1133 if (is_mandatory(i) .and. .not. got_this_field(i)) then
1134 call mprintf(.true.,ERROR,'The mandatory field %s was not found in any input data.',s1=fieldname(i))
1135 end if
1136 end do
1138 !
1139 ! Before we begin to write fields, if debug_level is set high enough, we
1140 ! write a table of which fields are available at which levels to the
1141 ! metgrid.log file, and then we check to see if any fields are not
1142 ! completely covered with data.
1143 !
1144 call storage_print_fields()
1145 call find_missing_values()
1147 !
1148 ! All of the processing is now done for this time period for this domain;
1149 ! now we simply output every field from the storage module.
1150 !
1152 title = 'OUTPUT FROM METGRID V3.4.1'
1154 ! Initialize the output module for this domain and time
1155 call mprintf(.true.,LOGFILE,'Initializing output module.')
1156 output_date = temp_date
1157 if ( .not. nocolons ) then
1158 if (len_trim(temp_date) == 13) then
1159 output_date(14:19) = ':00:00'
1160 else if (len_trim(temp_date) == 16) then
1161 output_date(17:19) = ':00'
1162 end if
1163 else
1164 if (len_trim(temp_date) == 13) then
1165 output_date(14:19) = '_00_00'
1166 else if (len_trim(temp_date) == 16) then
1167 output_date(17:19) = '_00'
1168 end if
1169 endif
1171 call output_init(n, title, output_date, gridtype, dyn_opt, &
1172 corner_lats, corner_lons, &
1173 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1174 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, &
1175 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1176 extra_col, extra_row)
1178 call get_bottom_top_dim(bottom_top_dim)
1180 ! Add in a flag to tell real that we have seven new msf fields
1181 nmet_flags = num_entries + 1
1182 allocate(met_flags(nmet_flags))
1183 do i=1,num_entries
1184 met_flags(i) = output_flags(i)
1185 end do
1186 if (gridtype == 'C') then
1187 met_flags(num_entries+1) = 'FLAG_MF_XY'
1188 else
1189 met_flags(num_entries+1) = ' '
1190 end if
1192 ! Find out how many soil levels or layers we have; this assumes a field named ST
1193 field % header % field = 'ST'
1194 nullify(soil_levels)
1195 call storage_get_levels(field, soil_levels)
1197 if (.not. associated(soil_levels)) then
1198 field % header % field = 'SOILT'
1199 nullify(soil_levels)
1200 call storage_get_levels(field, soil_levels)
1201 end if
1203 if (.not. associated(soil_levels)) then
1204 field % header % field = 'STC_WPS'
1205 nullify(soil_levels)
1206 call storage_get_levels(field, soil_levels)
1207 end if
1209 if (associated(soil_levels)) then
1210 num_metgrid_soil_levs = size(soil_levels)
1211 deallocate(soil_levels)
1212 else
1213 num_metgrid_soil_levs = 0
1214 end if
1216 ! First write out global attributes
1217 call mprintf(.true.,LOGFILE,'Writing global attributes to output.')
1218 call write_global_attrs(title, output_date, gridtype, dyn_opt, west_east_dim, &
1219 south_north_dim, bottom_top_dim, &
1220 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
1221 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
1222 map_proj, mminlu, num_land_cat, &
1223 is_water, is_lake, is_ice, is_urban, i_soilwater, &
1224 grid_id, parent_id, i_parent_start, &
1225 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
1226 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
1227 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
1228 corner_lats, corner_lons, num_metgrid_soil_levs, &
1229 met_flags, nmet_flags, process_bdy_width)
1231 deallocate(met_flags)
1233 call reset_next_field()
1235 istatus = 0
1237 ! Now loop over all output fields, writing each to the output module
1238 do while (istatus == 0)
1239 call get_next_output_field(cname, real_array, &
1240 sm1, em1, sm2, em2, sm3, em3, istatus)
1241 if (istatus == 0) then
1243 call mprintf(.true.,LOGFILE,'Writing field %s to output.',s1=cname)
1244 call write_field(sm1, em1, sm2, em2, sm3, em3, &
1245 cname, output_date, real_array)
1246 deallocate(real_array)
1248 end if
1250 end do
1252 call mprintf(.true.,LOGFILE,'Closing output file.')
1253 call output_close()
1255 ! Free up memory used by met fields for this valid time
1256 call storage_delete_all_td()
1258 end subroutine process_single_met_time
1261 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1262 ! Name: get_interp_masks
1263 !
1264 ! Purpose:
1265 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1266 subroutine get_interp_masks(fg_prefix, is_constants, fg_date)
1268 use interp_option_module
1269 use read_met_module
1270 use storage_module
1272 implicit none
1274 ! Arguments
1275 logical, intent(in) :: is_constants
1276 character (len=*), intent(in) :: fg_prefix, fg_date
1278 ! BUG: Move this constant to misc_definitions_module?
1279 integer, parameter :: BDR_WIDTH = 3
1281 ! Local variables
1282 integer :: i, istatus, idx, idxt
1283 type (fg_input) :: mask_field
1284 type (met_data) :: fg_data
1286 istatus = 0
1288 call read_met_init(fg_prefix, is_constants, fg_date, istatus)
1290 do while (istatus == 0)
1292 call read_next_met_field(fg_data, istatus)
1294 if (istatus == 0) then
1296 ! Find out which METGRID.TBL entry goes with this field
1297 idxt = num_entries + 1
1298 do idx=1,num_entries
1299 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
1300 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
1302 if (index(fg_prefix,trim(from_input(idx))) /= 0 .or. &
1303 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1304 idxt = idx
1305 end if
1307 end if
1308 end do
1309 idx = idxt
1310 if (idx > num_entries) idx = num_entries ! The last entry is a default
1312 ! Do we need to rename this field?
1313 if (output_name(idx) /= ' ') then
1314 fg_data%field = output_name(idx)(1:9)
1316 idxt = num_entries + 1
1317 do idx=1,num_entries
1318 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
1319 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
1321 if (index(fg_prefix,trim(from_input(idx))) /= 0 .or. &
1322 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1323 idxt = idx
1324 end if
1326 end if
1327 end do
1328 idx = idxt
1329 if (idx > num_entries) idx = num_entries ! The last entry is a default
1330 end if
1332 do i=1,num_entries
1333 if (interp_mask(i) /= ' ' .and. (trim(interp_mask(i)) == trim(fg_data%field))) then
1335 mask_field%header%version = 1
1336 mask_field%header%date = ' '
1337 mask_field%header%date = fg_date
1338 if (is_constants) then
1339 mask_field%header%time_dependent = .false.
1340 else
1341 mask_field%header%time_dependent = .true.
1342 end if
1343 mask_field%header%mask_field = .true.
1344 mask_field%header%forecast_hour = 0.
1345 mask_field%header%fg_source = 'degribbed met data'
1346 mask_field%header%field = trim(fg_data%field)//'.mask'
1347 mask_field%header%units = '-'
1348 mask_field%header%description = '-'
1349 mask_field%header%vertical_coord = 'none'
1350 mask_field%header%vertical_level = 1
1351 mask_field%header%sr_x = 1
1352 mask_field%header%sr_y = 1
1353 mask_field%header%array_order = 'XY'
1354 mask_field%header%dim1(1) = 1
1355 mask_field%header%dim1(2) = fg_data%nx
1356 mask_field%header%dim2(1) = 1
1357 mask_field%header%dim2(2) = fg_data%ny
1358 mask_field%header%is_wind_grid_rel = .true.
1359 mask_field%header%array_has_missing_values = .false.
1360 mask_field%map%stagger = M
1362 ! Do a simple check to see whether this is a global lat/lon dataset
1363 ! Note that we do not currently support regional Gaussian grids
1364 if ((fg_data%iproj == PROJ_LATLON .and. abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) &
1365 .or. (fg_data%iproj == PROJ_GAUSS)) then
1366 allocate(mask_field%r_arr(1-BDR_WIDTH:fg_data%nx+BDR_WIDTH,1:fg_data%ny))
1368 mask_field%r_arr(1:fg_data%nx, 1:fg_data%ny) = &
1369 fg_data%slab(1:fg_data%nx, 1:fg_data%ny)
1371 mask_field%r_arr(1-BDR_WIDTH:0, 1:fg_data%ny) = &
1372 fg_data%slab(fg_data%nx-BDR_WIDTH+1:fg_data%nx, 1:fg_data%ny)
1374 mask_field%r_arr(fg_data%nx+1:fg_data%nx+BDR_WIDTH, 1:fg_data%ny) = &
1375 fg_data%slab(1:BDR_WIDTH, 1:fg_data%ny)
1377 else
1378 allocate(mask_field%r_arr(1:fg_data%nx,1:fg_data%ny))
1379 mask_field%r_arr = fg_data%slab
1380 end if
1382 nullify(mask_field%valid_mask)
1383 nullify(mask_field%modified_mask)
1385 call storage_put_field(mask_field)
1387 exit
1389 end if
1390 end do
1392 if (associated(fg_data%slab)) deallocate(fg_data%slab)
1394 end if
1396 end do
1398 call read_met_close()
1400 end subroutine get_interp_masks
1403 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1404 ! Name: interp_met_field
1405 !
1406 ! Purpose:
1407 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1408 subroutine interp_met_field(input_name, short_fieldnm, ifieldstagger, istagger, &
1409 field, xlat, xlon, sm1, em1, sm2, em2, &
1410 sd1, ed1, sd2, ed2, &
1411 slab, minx, maxx, miny, maxy, bdr, do_gcell_interp, &
1412 new_pts, process_bdy_width, landmask)
1414 use bitarray_module
1415 use interp_module
1416 use interp_option_module
1417 use llxy_module
1418 use misc_definitions_module
1419 use storage_module
1421 implicit none
1423 ! Arguments
1424 integer, intent(in) :: ifieldstagger, istagger, &
1425 sm1, em1, sm2, em2, &
1426 sd1, ed1, sd2, ed2, &
1427 minx, maxx, miny, maxy, bdr, &
1428 process_bdy_width
1429 real, dimension(minx:maxx,miny:maxy), intent(in) :: slab
1430 real, dimension(sm1:em1,sm2:em2), intent(in) :: xlat, xlon
1431 real, dimension(sm1:em1,sm2:em2), intent(in), optional :: landmask
1432 logical, intent(in) :: do_gcell_interp
1433 character (len=9), intent(in) :: short_fieldnm
1434 character (len=MAX_FILENAME_LEN), intent(in) :: input_name
1435 type (fg_input), intent(inout) :: field
1436 type (bitarray), intent(inout) :: new_pts
1438 ! Local variables
1439 integer :: i, j, idx, idxt, orig_selected_proj, interp_mask_status, &
1440 interp_land_mask_status, interp_water_mask_status, process_width
1441 integer, pointer, dimension(:) :: interp_array
1442 real :: rx, ry, temp
1443 real, pointer, dimension(:,:) :: data_count
1444 type (fg_input) :: mask_field, mask_water_field, mask_land_field
1446 ! CWH Initialize local pointer variables
1447 nullify(interp_array)
1448 nullify(data_count)
1450 ! Find index into fieldname, interp_method, masked, and fill_missing
1451 ! of the current field
1452 idxt = num_entries + 1
1453 do idx=1,num_entries
1454 if ((index(fieldname(idx), trim(short_fieldnm)) /= 0) .and. &
1455 (len_trim(fieldname(idx)) == len_trim(short_fieldnm))) then
1456 if (index(input_name,trim(from_input(idx))) /= 0 .or. &
1457 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1458 idxt = idx
1459 end if
1460 end if
1461 end do
1462 idx = idxt
1463 if (idx > num_entries) then
1464 call mprintf(.true.,WARN,'Entry in METGRID.TBL not found for field %s. '// &
1465 'Default options will be used for this field!', s1=short_fieldnm)
1466 idx = num_entries ! The last entry is a default
1467 end if
1469 if (process_bdy_width == 0) then
1470 process_width = max(ed1-sd1+1, ed2-sd2+1)
1471 else
1472 process_width = process_bdy_width
1473 ! Add two extra rows/cols to accommodate staggered fields: one extra row/col for
1474 ! averaging to mass points in real, and one beyond that for averaging during
1475 ! wind rotation
1476 if (ifieldstagger /= M) process_width = process_width + 2
1477 end if
1479 field%header%dim1(1) = sm1
1480 field%header%dim1(2) = em1
1481 field%header%dim2(1) = sm2
1482 field%header%dim2(2) = em2
1483 field%map%stagger = ifieldstagger
1484 if (.not. associated(field%r_arr)) then
1485 allocate(field%r_arr(sm1:em1,sm2:em2))
1486 end if
1488 interp_mask_status = 1
1489 interp_land_mask_status = 1
1490 interp_water_mask_status = 1
1492 if (interp_mask(idx) /= ' ') then
1493 mask_field%header%version = 1
1494 mask_field%header%forecast_hour = 0.
1495 mask_field%header%field = trim(interp_mask(idx))//'.mask'
1496 mask_field%header%vertical_coord = 'none'
1497 mask_field%header%vertical_level = 1
1499 call storage_get_field(mask_field, interp_mask_status)
1501 end if
1502 if (interp_land_mask(idx) /= ' ') then
1503 mask_land_field%header%version = 1
1504 mask_land_field%header%forecast_hour = 0.
1505 mask_land_field%header%field = trim(interp_land_mask(idx))//'.mask'
1506 mask_land_field%header%vertical_coord = 'none'
1507 mask_land_field%header%vertical_level = 1
1509 call storage_get_field(mask_land_field, interp_land_mask_status)
1511 end if
1512 if (interp_water_mask(idx) /= ' ') then
1513 mask_water_field%header%version = 1
1514 mask_water_field%header%forecast_hour = 0.
1515 mask_water_field%header%field = trim(interp_water_mask(idx))//'.mask'
1516 mask_water_field%header%vertical_coord = 'none'
1517 mask_water_field%header%vertical_level = 1
1519 call storage_get_field(mask_water_field, interp_water_mask_status)
1521 end if
1523 interp_array => interp_array_from_string(interp_method(idx))
1526 !
1527 ! Interpolate using average_gcell interpolation method
1528 !
1529 if (do_gcell_interp) then
1530 allocate(data_count(sm1:em1,sm2:em2))
1531 data_count = 0.
1533 if (interp_mask_status == 0) then
1534 call accum_continuous(slab, &
1535 minx, maxx, miny, maxy, 1, 1, bdr, &
1536 field%r_arr, data_count, &
1537 sm1, em1, sm2, em2, 1, 1, &
1538 istagger, &
1539 new_pts, missing_value(idx), interp_mask_val(idx), interp_mask_relational(idx), mask_field%r_arr)
1540 else
1541 call accum_continuous(slab, &
1542 minx, maxx, miny, maxy, 1, 1, bdr, &
1543 field%r_arr, data_count, &
1544 sm1, em1, sm2, em2, 1, 1, &
1545 istagger, &
1546 new_pts, missing_value(idx), -1.) ! The -1 is the maskval, but since we
1547 ! we do not give an optional mask, no
1548 ! no need to worry about -1s in data
1549 end if
1551 orig_selected_proj = iget_selected_domain()
1552 call select_domain(SOURCE_PROJ)
1553 do j=sm2,em2
1554 do i=sm1,em1
1556 if (abs(i - sd1) >= process_width .and. (abs(i - ed1) >= process_width) .and. &
1557 abs(j - sd2) >= process_width .and. (abs(j - ed2) >= process_width)) then
1558 field%r_arr(i,j) = fill_missing(idx)
1559 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1560 cycle
1561 end if
1563 if (present(landmask)) then
1565 if (landmask(i,j) /= masked(idx)) then
1566 if (data_count(i,j) > 0.) then
1567 field%r_arr(i,j) = field%r_arr(i,j) / data_count(i,j)
1568 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1569 else
1571 if (interp_mask_status == 0) then
1572 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1573 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1574 mask_relational=interp_mask_relational(idx), &
1575 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1576 else
1577 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1578 minx, maxx, miny, maxy, bdr, missing_value(idx))
1579 end if
1581 if (temp /= missing_value(idx)) then
1582 field%r_arr(i,j) = temp
1583 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1584 end if
1586 end if
1587 else
1588 field%r_arr(i,j) = fill_missing(idx)
1589 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1590 end if
1592 if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
1593 .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
1594 field%r_arr(i,j) = fill_missing(idx)
1596 ! Assume that if missing fill value is other than default, then user has asked
1597 ! to fill in any missing values, and we can consider this point to have
1598 ! received a valid value
1599 if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1600 end if
1602 else
1604 if (data_count(i,j) > 0.) then
1605 field%r_arr(i,j) = field%r_arr(i,j) / data_count(i,j)
1606 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1607 else
1609 if (interp_mask_status == 0) then
1610 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1611 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1612 mask_relational=interp_mask_relational(idx), &
1613 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1614 else
1615 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1616 minx, maxx, miny, maxy, bdr, missing_value(idx))
1617 end if
1619 if (temp /= missing_value(idx)) then
1620 field%r_arr(i,j) = temp
1621 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1622 end if
1624 end if
1626 if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
1627 .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
1628 field%r_arr(i,j) = fill_missing(idx)
1630 ! Assume that if missing fill value is other than default, then user has asked
1631 ! to fill in any missing values, and we can consider this point to have
1632 ! received a valid value
1633 if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1634 end if
1636 end if
1638 end do
1639 end do
1640 call select_domain(orig_selected_proj)
1641 deallocate(data_count)
1643 !
1644 ! No average_gcell interpolation method
1645 !
1646 else
1648 orig_selected_proj = iget_selected_domain()
1649 call select_domain(SOURCE_PROJ)
1650 do j=sm2,em2
1651 do i=sm1,em1
1653 if (abs(i - sd1) >= process_width .and. (abs(i - ed1) >= process_width) .and. &
1654 abs(j - sd2) >= process_width .and. (abs(j - ed2) >= process_width)) then
1655 field%r_arr(i,j) = fill_missing(idx)
1656 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1657 cycle
1658 end if
1660 if (present(landmask)) then
1662 if (masked(idx) == MASKED_BOTH) then
1664 if (landmask(i,j) == 0) then ! WATER POINT
1666 if (interp_land_mask_status == 0) then
1667 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1668 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1669 mask_relational=interp_land_mask_relational(idx), &
1670 mask_val=interp_land_mask_val(idx), mask_field=mask_land_field%r_arr)
1671 else
1672 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1673 minx, maxx, miny, maxy, bdr, missing_value(idx))
1674 end if
1676 else if (landmask(i,j) == 1) then ! LAND POINT
1678 if (interp_water_mask_status == 0) then
1679 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1680 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1681 mask_relational=interp_water_mask_relational(idx), &
1682 mask_val=interp_water_mask_val(idx), mask_field=mask_water_field%r_arr)
1683 else
1684 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1685 minx, maxx, miny, maxy, bdr, missing_value(idx))
1686 end if
1688 end if
1690 else if (landmask(i,j) /= masked(idx)) then
1692 if (interp_mask_status == 0) then
1693 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, 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, slab, &
1699 minx, maxx, miny, maxy, bdr, missing_value(idx))
1700 end if
1702 else
1703 temp = missing_value(idx)
1704 end if
1706 ! No landmask for this field
1707 else
1709 if (interp_mask_status == 0) then
1710 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1711 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1712 mask_relational=interp_mask_relational(idx), &
1713 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1714 else
1715 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1716 minx, maxx, miny, maxy, bdr, missing_value(idx))
1717 end if
1719 end if
1721 if (temp /= missing_value(idx)) then
1722 field%r_arr(i,j) = temp
1723 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1724 else if (present(landmask)) then
1725 if (landmask(i,j) == masked(idx)) then
1726 field%r_arr(i,j) = fill_missing(idx)
1727 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1728 end if
1729 end if
1731 if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
1732 .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
1733 field%r_arr(i,j) = fill_missing(idx)
1735 ! Assume that if missing fill value is other than default, then user has asked
1736 ! to fill in any missing values, and we can consider this point to have
1737 ! received a valid value
1738 if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1739 end if
1741 end do
1742 end do
1743 call select_domain(orig_selected_proj)
1744 end if
1746 deallocate(interp_array)
1748 end subroutine interp_met_field
1751 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1752 ! Name: interp_to_latlon
1753 !
1754 ! Purpose:
1755 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1756 function interp_to_latlon(rlat, rlon, istagger, interp_method_list, slab, &
1757 minx, maxx, miny, maxy, bdr, source_missing_value, &
1758 mask_field, mask_relational, mask_val)
1760 use interp_module
1761 use llxy_module
1763 implicit none
1765 ! Arguments
1766 integer, intent(in) :: minx, maxx, miny, maxy, bdr, istagger
1767 integer, dimension(:), intent(in) :: interp_method_list
1768 real, intent(in) :: rlat, rlon, source_missing_value
1769 real, dimension(minx:maxx,miny:maxy), intent(in) :: slab
1770 real, intent(in), optional :: mask_val
1771 real, dimension(minx:maxx,miny:maxy), intent(in), optional :: mask_field
1772 character(len=1), intent(in), optional :: mask_relational
1774 ! Return value
1775 real :: interp_to_latlon
1777 ! Local variables
1778 real :: rx, ry
1780 interp_to_latlon = source_missing_value
1782 call lltoxy(rlat, rlon, rx, ry, istagger)
1783 if (rx >= minx+bdr-0.5 .and. rx <= maxx-bdr+0.5) then
1784 if (present(mask_field) .and. present(mask_val) .and. present (mask_relational)) then
1785 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
1786 interp_method_list, 1, mask_relational, mask_val, mask_field)
1787 else if (present(mask_field) .and. present(mask_val)) then
1788 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
1789 interp_method_list, 1, maskval=mask_val, mask_array=mask_field)
1790 else
1791 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
1792 interp_method_list, 1)
1793 end if
1794 else
1795 interp_to_latlon = source_missing_value
1796 end if
1798 if (interp_to_latlon == source_missing_value) then
1800 ! Try a lon in the range 0. to 360.; all lons in the xlon
1801 ! array should be in the range -180. to 180.
1802 if (rlon < 0.) then
1803 call lltoxy(rlat, rlon+360., rx, ry, istagger)
1804 if (rx >= minx+bdr-0.5 .and. rx <= maxx-bdr+0.5) then
1805 if (present(mask_field) .and. present(mask_val) .and. present(mask_relational)) then
1806 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
1807 1, 1, source_missing_value, &
1808 interp_method_list, 1, mask_relational, mask_val, mask_field)
1809 else if (present(mask_field) .and. present(mask_val)) then
1810 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
1811 1, 1, source_missing_value, &
1812 interp_method_list, 1, maskval=mask_val, mask_array=mask_field)
1813 else
1814 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
1815 1, 1, source_missing_value, &
1816 interp_method_list, 1)
1817 end if
1818 else
1819 interp_to_latlon = source_missing_value
1820 end if
1822 end if
1824 end if
1826 return
1828 end function interp_to_latlon
1831 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1832 ! Name: get_bottom_top_dim
1833 !
1834 ! Purpose:
1835 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1836 subroutine get_bottom_top_dim(bottom_top_dim)
1838 use interp_option_module
1839 use list_module
1840 use storage_module
1842 implicit none
1844 ! Arguments
1845 integer, intent(out) :: bottom_top_dim
1847 ! Local variables
1848 integer :: i, j
1849 integer, pointer, dimension(:) :: field_levels
1850 character (len=32) :: z_dim
1851 type (fg_input), pointer, dimension(:) :: headers
1852 type (list) :: temp_levels
1854 !CWH Initialize local pointer variables
1855 nullify(field_levels)
1856 nullify(headers)
1858 ! Initialize a list to store levels that are found for 3-d fields
1859 call list_init(temp_levels)
1861 ! Get a list of all time-dependent fields (given by their headers) from
1862 ! the storage module
1863 call storage_get_td_headers(headers)
1865 !
1866 ! Given headers of all fields, we first build a list of all possible levels
1867 ! for 3-d met fields (excluding sea-level, though).
1868 !
1869 do i=1,size(headers)
1870 call get_z_dim_name(headers(i)%header%field, z_dim)
1872 ! We only want to consider 3-d met fields
1873 if (z_dim(1:18) == 'num_metgrid_levels') then
1875 ! Find out what levels the current field has
1876 call storage_get_levels(headers(i), field_levels)
1877 do j=1,size(field_levels)
1879 ! If this level has not yet been encountered, add it to our list
1880 if (.not. list_search(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))) then
1881 if (field_levels(j) /= 201300) then
1882 call list_insert(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))
1883 end if
1884 end if
1886 end do
1888 deallocate(field_levels)
1890 end if
1892 end do
1894 bottom_top_dim = list_length(temp_levels)
1896 call list_destroy(temp_levels)
1897 deallocate(headers)
1899 end subroutine get_bottom_top_dim
1902 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1903 ! Name: fill_missing_levels
1904 !
1905 ! Purpose:
1906 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1907 subroutine fill_missing_levels(output_flags)
1909 use interp_option_module
1910 use list_module
1911 use module_debug
1912 use module_mergesort
1913 use storage_module
1915 implicit none
1917 ! Arguments
1918 character (len=128), dimension(:), intent(inout) :: output_flags
1920 ! Local variables
1921 integer :: i, ii, j, ix, jx, k, lower, upper, temp, istatus
1922 integer, pointer, dimension(:) :: union_levels, field_levels
1923 real, pointer, dimension(:) :: r_union_levels
1924 character (len=128) :: clevel
1925 type (fg_input) :: lower_field, upper_field, new_field, search_field
1926 type (fg_input), pointer, dimension(:) :: headers, all_headers
1927 type (list) :: temp_levels
1928 type (list_item), pointer, dimension(:) :: keys
1930 ! CWH Initialize local pointer variables
1931 nullify(union_levels)
1932 nullify(field_levels)
1933 nullify(r_union_levels)
1934 nullify(headers)
1935 nullify(all_headers)
1936 nullify(keys)
1938 ! Initialize a list to store levels that are found for 3-d fields
1939 call list_init(temp_levels)
1941 ! Get a list of all fields (given by their headers) from the storage module
1942 call storage_get_td_headers(headers)
1943 call storage_get_all_headers(all_headers)
1945 !
1946 ! Given headers of all fields, we first build a list of all possible levels
1947 ! for 3-d met fields (excluding sea-level, though).
1948 !
1949 do i=1,size(headers)
1951 ! Find out what levels the current field has
1952 call storage_get_levels(headers(i), field_levels)
1953 do j=1,size(field_levels)
1955 ! If this level has not yet been encountered, add it to our list
1956 if (.not. list_search(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))) then
1957 if (field_levels(j) /= 201300) then
1958 call list_insert(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))
1959 end if
1960 end if
1962 end do
1964 deallocate(field_levels)
1966 end do
1968 if (list_length(temp_levels) > 0) then
1970 !
1971 ! With all possible levels stored in a list, get an array of levels, sorted
1972 ! in decreasing order
1973 !
1974 i = 0
1975 allocate(union_levels(list_length(temp_levels)))
1976 do while (list_length(temp_levels) > 0)
1977 i = i + 1
1978 call list_get_first_item(temp_levels, ikey=union_levels(i), ivalue=temp)
1979 end do
1980 call mergesort(union_levels, 1, size(union_levels))
1982 allocate(r_union_levels(size(union_levels)))
1983 do i=1,size(union_levels)
1984 r_union_levels(i) = real(union_levels(i))
1985 end do
1987 !
1988 ! With a sorted, complete list of levels, we need
1989 ! to go back and fill in missing levels for each 3-d field
1990 !
1991 do i=1,size(headers)
1993 !
1994 ! Find entry in METGRID.TBL for this field, if one exists; if it does, then the
1995 ! entry may tell us how to get values for the current field at the missing level
1996 !
1997 do ii=1,num_entries
1998 if (fieldname(ii) == headers(i)%header%field) exit
1999 end do
2000 if (ii <= num_entries) then
2001 call dup(headers(i),new_field)
2002 nullify(new_field%valid_mask)
2003 nullify(new_field%modified_mask)
2004 call fill_field(new_field, ii, output_flags, r_union_levels)
2005 end if
2007 end do
2009 deallocate(union_levels)
2010 deallocate(r_union_levels)
2011 deallocate(headers)
2013 call storage_get_td_headers(headers)
2015 !
2016 ! Now we may need to vertically interpolate to missing values in 3-d fields
2017 !
2018 do i=1,size(headers)
2020 call storage_get_levels(headers(i), field_levels)
2022 ! If this isn't a 3-d array, nothing to do
2023 if (size(field_levels) > 1) then
2025 do k=1,size(field_levels)
2026 call dup(headers(i),search_field)
2027 search_field%header%vertical_level = field_levels(k)
2028 call storage_get_field(search_field,istatus)
2029 if (istatus == 0) then
2030 JLOOP: do jx=search_field%header%dim2(1),search_field%header%dim2(2)
2031 ILOOP: do ix=search_field%header%dim1(1),search_field%header%dim1(2)
2032 if (.not. bitarray_test(search_field%valid_mask, &
2033 ix-search_field%header%dim1(1)+1, &
2034 jx-search_field%header%dim2(1)+1)) then
2036 call dup(search_field, lower_field)
2037 do lower=k-1,1,-1
2038 lower_field%header%vertical_level = field_levels(lower)
2039 call storage_get_field(lower_field,istatus)
2040 if (bitarray_test(lower_field%valid_mask, &
2041 ix-search_field%header%dim1(1)+1, &
2042 jx-search_field%header%dim2(1)+1)) &
2043 exit
2045 end do
2047 call dup(search_field, upper_field)
2048 do upper=k+1,size(field_levels)
2049 upper_field%header%vertical_level = field_levels(upper)
2050 call storage_get_field(upper_field,istatus)
2051 if (bitarray_test(upper_field%valid_mask, &
2052 ix-search_field%header%dim1(1)+1, &
2053 jx-search_field%header%dim2(1)+1)) &
2054 exit
2056 end do
2057 if (upper <= size(field_levels) .and. lower >= 1) then
2058 search_field%r_arr(ix,jx) = real(abs(field_levels(upper)-field_levels(k))) &
2059 / real(abs(field_levels(upper)-field_levels(lower))) &
2060 * lower_field%r_arr(ix,jx) &
2061 + real(abs(field_levels(k)-field_levels(lower))) &
2062 / real(abs(field_levels(upper)-field_levels(lower))) &
2063 * upper_field%r_arr(ix,jx)
2064 call bitarray_set(search_field%valid_mask, &
2065 ix-search_field%header%dim1(1)+1, &
2066 jx-search_field%header%dim2(1)+1)
2067 end if
2068 end if
2069 end do ILOOP
2070 end do JLOOP
2071 else
2072 call mprintf(.true.,ERROR, &
2073 'This is bad, could not get %s at level %i.', &
2074 s1=trim(search_field%header%field), i1=field_levels(k))
2075 end if
2076 end do
2078 end if
2080 deallocate(field_levels)
2082 end do
2084 end if
2086 call list_destroy(temp_levels)
2087 deallocate(all_headers)
2088 deallocate(headers)
2090 end subroutine fill_missing_levels
2093 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2094 ! Name: create_derived_fields
2095 !
2096 ! Purpose:
2097 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2098 subroutine create_derived_fields(arg_gridtype, hdate, xfcst, &
2099 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
2100 we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e, &
2101 created_this_field, output_flags)
2103 use interp_option_module
2104 use list_module
2105 use module_mergesort
2106 use storage_module
2108 implicit none
2110 ! Arguments
2111 integer, intent(in) :: we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
2112 we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e
2113 real, intent(in) :: xfcst
2114 logical, dimension(:), intent(inout) :: created_this_field
2115 character (len=1), intent(in) :: arg_gridtype
2116 character (len=24), intent(in) :: hdate
2117 character (len=128), dimension(:), intent(inout) :: output_flags
2119 ! Local variables
2120 integer :: idx, i, j, istatus
2121 type (fg_input) :: field
2123 ! Initialize fg_input structure to store the field
2124 field%header%version = 1
2125 field%header%date = hdate//' '
2126 field%header%time_dependent = .true.
2127 field%header%mask_field = .false.
2128 field%header%forecast_hour = xfcst
2129 field%header%fg_source = 'Derived from FG'
2130 field%header%field = ' '
2131 field%header%units = ' '
2132 field%header%description = ' '
2133 field%header%vertical_level = 0
2134 field%header%sr_x = 1
2135 field%header%sr_y = 1
2136 field%header%array_order = 'XY '
2137 field%header%is_wind_grid_rel = .true.
2138 field%header%array_has_missing_values = .false.
2139 nullify(field%r_arr)
2140 nullify(field%valid_mask)
2141 nullify(field%modified_mask)
2143 !
2144 ! Check each entry in METGRID.TBL to see whether it is a derive field
2145 !
2146 do idx=1,num_entries
2147 if (is_derived_field(idx)) then
2149 created_this_field(idx) = .true.
2151 call mprintf(.true.,INFORM,'Going to create the field %s',s1=fieldname(idx))
2153 ! Intialize more fields in storage structure
2154 field%header%field = fieldname(idx)
2155 call get_z_dim_name(fieldname(idx),field%header%vertical_coord)
2156 field%map%stagger = output_stagger(idx)
2157 if (arg_gridtype == 'E') then
2158 field%header%dim1(1) = we_mem_s
2159 field%header%dim1(2) = we_mem_e
2160 field%header%dim2(1) = sn_mem_s
2161 field%header%dim2(2) = sn_mem_e
2162 else if (arg_gridtype == 'C') then
2163 if (output_stagger(idx) == M) then
2164 field%header%dim1(1) = we_mem_s
2165 field%header%dim1(2) = we_mem_e
2166 field%header%dim2(1) = sn_mem_s
2167 field%header%dim2(2) = sn_mem_e
2168 else if (output_stagger(idx) == U) then
2169 field%header%dim1(1) = we_mem_stag_s
2170 field%header%dim1(2) = we_mem_stag_e
2171 field%header%dim2(1) = sn_mem_s
2172 field%header%dim2(2) = sn_mem_e
2173 else if (output_stagger(idx) == V) then
2174 field%header%dim1(1) = we_mem_s
2175 field%header%dim1(2) = we_mem_e
2176 field%header%dim2(1) = sn_mem_stag_s
2177 field%header%dim2(2) = sn_mem_stag_e
2178 end if
2179 end if
2181 call fill_field(field, idx, output_flags)
2183 end if
2184 end do
2187 end subroutine create_derived_fields
2190 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2191 ! Name: fill_field
2192 !
2193 ! Purpose:
2194 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2195 subroutine fill_field(field, idx, output_flags, all_level_list)
2197 use interp_option_module
2198 use list_module
2199 use module_mergesort
2200 use storage_module
2202 implicit none
2204 ! Arguments
2205 integer, intent(in) :: idx
2206 type (fg_input), intent(inout) :: field
2207 character (len=128), dimension(:), intent(inout) :: output_flags
2208 real, dimension(:), intent(in), optional :: all_level_list
2210 ! Local variables
2211 integer :: i, j, istatus, isrclevel
2212 integer, pointer, dimension(:) :: all_list
2213 real :: rfillconst, rlevel, rsrclevel
2214 type (fg_input) :: query_field
2215 type (list_item), pointer, dimension(:) :: keys
2216 character (len=128) :: asrcname
2217 logical :: filled_all_lev
2219 !CWH Initialize local pointer variables
2220 nullify(all_list)
2221 nullify(keys)
2223 filled_all_lev = .false.
2225 !
2226 ! Get a list of all levels to be filled for this field
2227 !
2228 keys => list_get_keys(fill_lev_list(idx))
2230 do i=1,list_length(fill_lev_list(idx))
2232 !
2233 ! First handle a specification for levels "all"
2234 !
2235 if (trim(keys(i)%ckey) == 'all') then
2237 ! We only want to fill all levels if we haven't already filled "all" of them
2238 if (.not. filled_all_lev) then
2240 filled_all_lev = .true.
2242 query_field%header%time_dependent = .true.
2243 query_field%header%field = ' '
2244 nullify(query_field%r_arr)
2245 nullify(query_field%valid_mask)
2246 nullify(query_field%modified_mask)
2248 ! See if we are filling this level with a constant
2249 call get_constant_fill_lev(keys(i)%cvalue, rfillconst, istatus)
2250 if (istatus == 0) then
2251 if (present(all_level_list)) then
2252 do j=1,size(all_level_list)
2253 call create_level(field, real(all_level_list(j)), idx, output_flags, rfillconst=rfillconst)
2254 end do
2255 else
2256 query_field%header%field = level_template(idx)
2257 nullify(all_list)
2258 call storage_get_levels(query_field, all_list)
2259 if (associated(all_list)) then
2260 do j=1,size(all_list)
2261 call create_level(field, real(all_list(j)), idx, output_flags, rfillconst=rfillconst)
2262 end do
2263 deallocate(all_list)
2264 end if
2265 end if
2267 ! Else see if we are filling this level with a constant equal
2268 ! to the value of the level
2269 else if (trim(keys(i)%cvalue) == 'vertical_index') then
2270 if (present(all_level_list)) then
2271 do j=1,size(all_level_list)
2272 call create_level(field, real(all_level_list(j)), idx, output_flags, &
2273 rfillconst=real(all_level_list(j)))
2274 end do
2275 else
2276 query_field%header%field = level_template(idx)
2277 nullify(all_list)
2278 call storage_get_levels(query_field, all_list)
2279 if (associated(all_list)) then
2280 do j=1,size(all_list)
2281 call create_level(field, real(all_list(j)), idx, output_flags, rfillconst=real(all_list(j)))
2282 end do
2283 deallocate(all_list)
2284 end if
2285 end if
2287 ! Else, we assume that it is a field from which we are copying levels
2288 else
2289 if (present(all_level_list)) then
2290 do j=1,size(all_level_list)
2291 call create_level(field, real(all_level_list(j)), idx, output_flags, &
2292 asrcname=keys(i)%cvalue, rsrclevel=real(all_level_list(j)))
2293 end do
2294 else
2295 query_field%header%field = keys(i)%cvalue ! Use same levels as source field, not level_template
2296 nullify(all_list)
2297 call storage_get_levels(query_field, all_list)
2298 if (associated(all_list)) then
2299 do j=1,size(all_list)
2300 call create_level(field, real(all_list(j)), idx, output_flags, &
2301 asrcname=keys(i)%cvalue, rsrclevel=real(all_list(j)))
2302 end do
2303 deallocate(all_list)
2305 else
2307 ! If the field doesn't have any levels (or does not exist) then we have not
2308 ! really filled all levels at this point.
2309 filled_all_lev = .false.
2310 end if
2311 end if
2313 end if
2314 end if
2316 !
2317 ! Handle individually specified levels
2318 !
2319 else
2321 read(keys(i)%ckey,*) rlevel
2323 ! See if we are filling this level with a constant
2324 call get_constant_fill_lev(keys(i)%cvalue, rfillconst, istatus)
2325 if (istatus == 0) then
2326 call create_level(field, rlevel, idx, output_flags, rfillconst=rfillconst)
2328 ! Otherwise, we are filling from another level
2329 else
2330 call get_fill_src_level(keys(i)%cvalue, asrcname, isrclevel)
2331 rsrclevel = real(isrclevel)
2332 call create_level(field, rlevel, idx, output_flags, &
2333 asrcname=asrcname, rsrclevel=rsrclevel)
2335 end if
2336 end if
2337 end do
2339 if (associated(keys)) deallocate(keys)
2341 end subroutine fill_field
2344 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2345 ! Name: create_level
2346 !
2347 ! Purpose:
2348 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2349 subroutine create_level(field_template, rlevel, idx, output_flags, &
2350 rfillconst, asrcname, rsrclevel)
2352 use storage_module
2353 use interp_option_module
2355 implicit none
2357 ! Arguments
2358 type (fg_input), intent(inout) :: field_template
2359 real, intent(in) :: rlevel
2360 integer, intent(in) :: idx
2361 character (len=128), dimension(:), intent(inout) :: output_flags
2362 real, intent(in), optional :: rfillconst, rsrclevel
2363 character (len=128), intent(in), optional :: asrcname
2365 ! Local variables
2366 integer :: i, j, istatus
2367 integer :: sm1,em1,sm2,em2
2368 type (fg_input) :: query_field
2370 !
2371 ! Check to make sure optional arguments are sane
2372 !
2373 if (present(rfillconst) .and. (present(asrcname) .or. present(rsrclevel))) then
2374 call mprintf(.true.,ERROR,'A call to create_level() cannot be given specifications '// &
2375 'for both a constant fill value and a source level.')
2377 else if ((present(asrcname) .and. .not. present(rsrclevel)) .or. &
2378 (.not. present(asrcname) .and. present(rsrclevel))) then
2379 call mprintf(.true.,ERROR,'Neither or both of optional arguments asrcname and '// &
2380 'rsrclevel must be specified to subroutine create_level().')
2382 else if (.not. present(rfillconst) .and. &
2383 .not. present(asrcname) .and. &
2384 .not. present(rsrclevel)) then
2385 call mprintf(.true.,ERROR,'A call to create_level() must be given either a specification '// &
2386 'for a constant fill value or a source level.')
2387 end if
2389 query_field%header%time_dependent = .true.
2390 query_field%header%field = field_template%header%field
2391 query_field%header%vertical_level = rlevel
2392 nullify(query_field%r_arr)
2393 nullify(query_field%valid_mask)
2394 nullify(query_field%modified_mask)
2396 call storage_query_field(query_field, istatus)
2397 if (istatus == 0) then
2398 call mprintf(.true.,INFORM,'%s at level %f already exists; leaving it alone.', &
2399 s1=field_template%header%field, f1=rlevel)
2400 return
2401 end if
2403 sm1 = field_template%header%dim1(1)
2404 em1 = field_template%header%dim1(2)
2405 sm2 = field_template%header%dim2(1)
2406 em2 = field_template%header%dim2(2)
2408 !
2409 ! Handle constant fill value case
2410 !
2411 if (present(rfillconst)) then
2413 field_template%header%vertical_level = rlevel
2414 allocate(field_template%r_arr(sm1:em1,sm2:em2))
2415 allocate(field_template%valid_mask)
2416 allocate(field_template%modified_mask)
2417 call bitarray_create(field_template%valid_mask, em1-sm1+1, em2-sm2+1)
2418 call bitarray_create(field_template%modified_mask, em1-sm1+1, em2-sm2+1)
2420 field_template%r_arr = rfillconst
2422 do j=sm2,em2
2423 do i=sm1,em1
2424 call bitarray_set(field_template%valid_mask, i-sm1+1, j-sm2+1)
2425 end do
2426 end do
2428 call storage_put_field(field_template)
2430 if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
2431 output_flags(idx) = flag_in_output(idx)
2432 end if
2434 !
2435 ! Handle source field and source level case
2436 !
2437 else if (present(asrcname) .and. present(rsrclevel)) then
2439 query_field%header%field = ' '
2440 query_field%header%field = asrcname
2441 query_field%header%vertical_level = rsrclevel
2443 ! Check to see whether the requested source field exists at the requested level
2444 call storage_query_field(query_field, istatus)
2446 if (istatus == 0) then
2448 ! Read in requested field at requested level
2449 call storage_get_field(query_field, istatus)
2450 if ((query_field%header%dim1(1) /= field_template%header%dim1(1)) .or. &
2451 (query_field%header%dim1(2) /= field_template%header%dim1(2)) .or. &
2452 (query_field%header%dim2(1) /= field_template%header%dim2(1)) .or. &
2453 (query_field%header%dim2(2) /= field_template%header%dim2(2))) then
2454 call mprintf(.true.,ERROR,'Dimensions for %s do not match those of %s. This is '// &
2455 'probably because the staggerings of the fields do not match.', &
2456 s1=query_field%header%field, s2=field_template%header%field)
2457 end if
2459 field_template%header%vertical_level = rlevel
2460 allocate(field_template%r_arr(sm1:em1,sm2:em2))
2461 allocate(field_template%valid_mask)
2462 allocate(field_template%modified_mask)
2463 call bitarray_create(field_template%valid_mask, em1-sm1+1, em2-sm2+1)
2464 call bitarray_create(field_template%modified_mask, em1-sm1+1, em2-sm2+1)
2466 field_template%r_arr = query_field%r_arr
2468 ! We should retain information about which points in the field are valid
2469 do j=sm2,em2
2470 do i=sm1,em1
2471 if (bitarray_test(query_field%valid_mask, i-sm1+1, j-sm2+1)) then
2472 call bitarray_set(field_template%valid_mask, i-sm1+1, j-sm2+1)
2473 end if
2474 end do
2475 end do
2477 call storage_put_field(field_template)
2479 if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
2480 output_flags(idx) = flag_in_output(idx)
2481 end if
2483 else
2484 call mprintf(.true.,INFORM,'Couldn''t find %s at level %f to fill level %f of %s.', &
2485 s1=asrcname,f1=rsrclevel,f2=rlevel,s2=field_template%header%field)
2486 end if
2488 end if
2490 end subroutine create_level
2493 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2494 ! Name: accum_continuous
2495 !
2496 ! Purpose: Sum up all of the source data points whose nearest neighbor in the
2497 ! model grid is the specified model grid point.
2498 !
2499 ! NOTE: When processing the source tile, those source points that are
2500 ! closest to a different model grid point will be added to the totals for
2501 ! such grid points; thus, an entire source tile will be processed at a time.
2502 ! This routine really processes for all model grid points that are
2503 ! within a source tile, and not just for a single grid point.
2504 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2505 subroutine accum_continuous(src_array, &
2506 src_min_x, src_max_x, src_min_y, src_max_y, src_min_z, src_max_z, bdr_width, &
2507 dst_array, n, &
2508 start_i, end_i, start_j, end_j, start_k, end_k, &
2509 istagger, &
2510 new_pts, msgval, maskval, mask_relational, mask_array, sr_x, sr_y)
2512 use bitarray_module
2513 use misc_definitions_module
2515 implicit none
2517 ! Arguments
2518 integer, intent(in) :: start_i, end_i, start_j, end_j, start_k, end_k, istagger, &
2519 src_min_x, src_max_x, src_min_y, src_max_y, src_min_z, src_max_z, bdr_width
2520 real, intent(in) :: maskval, msgval
2521 real, dimension(src_min_x:src_max_x, src_min_y:src_max_y, src_min_z:src_max_z), intent(in) :: src_array
2522 real, dimension(start_i:end_i, start_j:end_j, start_k:end_k), intent(inout) :: dst_array, n
2523 real, dimension(src_min_x:src_max_x, src_min_y:src_max_y), intent(in), optional :: mask_array
2524 integer, intent(in), optional :: sr_x, sr_y
2525 type (bitarray), intent(inout) :: new_pts
2526 character(len=1), intent(in), optional :: mask_relational
2528 ! Local variables
2529 integer :: i, j
2530 integer, pointer, dimension(:,:,:) :: where_maps_to
2531 real :: rsr_x, rsr_y
2533 rsr_x = 1.0
2534 rsr_y = 1.0
2535 if (present(sr_x)) rsr_x = real(sr_x)
2536 if (present(sr_y)) rsr_y = real(sr_y)
2538 allocate(where_maps_to(src_min_x:src_max_x,src_min_y:src_max_y,2))
2539 do i=src_min_x,src_max_x
2540 do j=src_min_y,src_max_y
2541 where_maps_to(i,j,1) = NOT_PROCESSED
2542 end do
2543 end do
2545 call process_continuous_block(src_array, where_maps_to, &
2546 src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
2547 src_min_x+bdr_width, src_min_y, src_min_z, &
2548 src_max_x-bdr_width, src_max_y, src_max_z, &
2549 dst_array, n, start_i, end_i, start_j, end_j, start_k, end_k, &
2550 istagger, &
2551 new_pts, rsr_x, rsr_y, msgval, maskval, mask_relational, mask_array)
2553 deallocate(where_maps_to)
2555 end subroutine accum_continuous
2558 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2559 ! Name: process_continuous_block
2560 !
2561 ! Purpose: To recursively process a subarray of continuous data, adding the
2562 ! points in a block to the sum for their nearest grid point. The nearest
2563 ! neighbor may be estimated in some cases; for example, if the four corners
2564 ! of a subarray all have the same nearest grid point, all elements in the
2565 ! subarray are added to that grid point.
2566 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2567 recursive subroutine process_continuous_block(tile_array, where_maps_to, &
2568 src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
2569 min_i, min_j, min_k, max_i, max_j, max_k, &
2570 dst_array, n, &
2571 start_x, end_x, start_y, end_y, start_z, end_z, &
2572 istagger, &
2573 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2575 use bitarray_module
2576 use llxy_module
2577 use misc_definitions_module
2579 implicit none
2581 ! Arguments
2582 integer, intent(in) :: min_i, min_j, min_k, max_i, max_j, max_k, &
2583 src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
2584 start_x, end_x, start_y, end_y, start_z, end_z, istagger
2585 integer, dimension(src_min_x:src_max_x,src_min_y:src_max_y,2), intent(inout) :: where_maps_to
2586 real, intent(in) :: sr_x, sr_y, maskval, msgval
2587 real, dimension(src_min_x:src_max_x,src_min_y:src_max_y,src_min_z:src_max_z), intent(in) :: tile_array
2588 real, dimension(src_min_x:src_max_x,src_min_y:src_max_y), intent(in), optional :: mask_array
2589 real, dimension(start_x:end_x,start_y:end_y,start_z:end_z), intent(inout) :: dst_array, n
2590 type (bitarray), intent(inout) :: new_pts
2591 character(len=1), intent(in), optional :: mask_relational
2593 ! Local variables
2594 integer :: orig_selected_domain, x_dest, y_dest, i, j, k, center_i, center_j
2595 real :: lat_corner, lon_corner, rx, ry
2597 ! Compute the model grid point that the corners of the rectangle to be
2598 ! processed map to
2599 ! Lower-left corner
2600 if (where_maps_to(min_i,min_j,1) == NOT_PROCESSED) then
2601 orig_selected_domain = iget_selected_domain()
2602 call select_domain(SOURCE_PROJ)
2603 call xytoll(real(min_i), real(min_j), lat_corner, lon_corner, istagger)
2604 call select_domain(1)
2605 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2606 rx = (rx - 1.0)*sr_x + 1.0
2607 ry = (ry - 1.0)*sr_y + 1.0
2608 call select_domain(orig_selected_domain)
2609 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2610 real(start_y) <= ry .and. ry <= real(end_y)) then
2611 where_maps_to(min_i,min_j,1) = nint(rx)
2612 where_maps_to(min_i,min_j,2) = nint(ry)
2613 else
2614 where_maps_to(min_i,min_j,1) = OUTSIDE_DOMAIN
2615 end if
2616 end if
2618 ! Upper-left corner
2619 if (where_maps_to(min_i,max_j,1) == NOT_PROCESSED) then
2620 orig_selected_domain = iget_selected_domain()
2621 call select_domain(SOURCE_PROJ)
2622 call xytoll(real(min_i), real(max_j), lat_corner, lon_corner, istagger)
2623 call select_domain(1)
2624 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2625 rx = (rx - 1.0)*sr_x + 1.0
2626 ry = (ry - 1.0)*sr_y + 1.0
2627 call select_domain(orig_selected_domain)
2628 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2629 real(start_y) <= ry .and. ry <= real(end_y)) then
2630 where_maps_to(min_i,max_j,1) = nint(rx)
2631 where_maps_to(min_i,max_j,2) = nint(ry)
2632 else
2633 where_maps_to(min_i,max_j,1) = OUTSIDE_DOMAIN
2634 end if
2635 end if
2637 ! Upper-right corner
2638 if (where_maps_to(max_i,max_j,1) == NOT_PROCESSED) then
2639 orig_selected_domain = iget_selected_domain()
2640 call select_domain(SOURCE_PROJ)
2641 call xytoll(real(max_i), real(max_j), lat_corner, lon_corner, istagger)
2642 call select_domain(1)
2643 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2644 rx = (rx - 1.0)*sr_x + 1.0
2645 ry = (ry - 1.0)*sr_y + 1.0
2646 call select_domain(orig_selected_domain)
2647 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2648 real(start_y) <= ry .and. ry <= real(end_y)) then
2649 where_maps_to(max_i,max_j,1) = nint(rx)
2650 where_maps_to(max_i,max_j,2) = nint(ry)
2651 else
2652 where_maps_to(max_i,max_j,1) = OUTSIDE_DOMAIN
2653 end if
2654 end if
2656 ! Lower-right corner
2657 if (where_maps_to(max_i,min_j,1) == NOT_PROCESSED) then
2658 orig_selected_domain = iget_selected_domain()
2659 call select_domain(SOURCE_PROJ)
2660 call xytoll(real(max_i), real(min_j), lat_corner, lon_corner, istagger)
2661 call select_domain(1)
2662 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2663 rx = (rx - 1.0)*sr_x + 1.0
2664 ry = (ry - 1.0)*sr_y + 1.0
2665 call select_domain(orig_selected_domain)
2666 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2667 real(start_y) <= ry .and. ry <= real(end_y)) then
2668 where_maps_to(max_i,min_j,1) = nint(rx)
2669 where_maps_to(max_i,min_j,2) = nint(ry)
2670 else
2671 where_maps_to(max_i,min_j,1) = OUTSIDE_DOMAIN
2672 end if
2673 end if
2675 ! If all four corners map to same model grid point, accumulate the
2676 ! entire rectangle
2677 if (where_maps_to(min_i,min_j,1) == where_maps_to(min_i,max_j,1) .and. &
2678 where_maps_to(min_i,min_j,1) == where_maps_to(max_i,max_j,1) .and. &
2679 where_maps_to(min_i,min_j,1) == where_maps_to(max_i,min_j,1) .and. &
2680 where_maps_to(min_i,min_j,2) == where_maps_to(min_i,max_j,2) .and. &
2681 where_maps_to(min_i,min_j,2) == where_maps_to(max_i,max_j,2) .and. &
2682 where_maps_to(min_i,min_j,2) == where_maps_to(max_i,min_j,2) .and. &
2683 where_maps_to(min_i,min_j,1) /= OUTSIDE_DOMAIN) then
2684 x_dest = where_maps_to(min_i,min_j,1)
2685 y_dest = where_maps_to(min_i,min_j,2)
2687 ! If this grid point was already given a value from higher-priority source data,
2688 ! there is nothing to do.
2689 ! if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2691 ! If this grid point has never been given a value by this level of source data,
2692 ! initialize the point
2693 if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2694 do k=min_k,max_k
2695 dst_array(x_dest,y_dest,k) = 0.
2696 end do
2697 end if
2699 ! Sum all the points whose nearest neighbor is this grid point
2700 if (present(mask_array) .and. present(mask_relational)) then
2701 do i=min_i,max_i
2702 do j=min_j,max_j
2703 do k=min_k,max_k
2704 ! Ignore masked/missing values in the source data
2705 if (tile_array(i,j,k) /= msgval) then
2706 if (mask_relational == ' ' .and. mask_array(i,j) /= maskval) then
2707 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2708 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2709 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2710 else if (mask_relational == '<' .and. mask_array(i,j) >= maskval) then
2711 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2712 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2713 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2714 else if (mask_relational == '>' .and. mask_array(i,j) <= maskval) then
2715 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2716 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2717 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2718 end if
2719 end if
2720 end do
2721 end do
2722 end do
2723 else if (present(mask_array)) then
2724 do i=min_i,max_i
2725 do j=min_j,max_j
2726 do k=min_k,max_k
2727 ! Ignore masked/missing values in the source data
2728 if ((tile_array(i,j,k) /= msgval) .and. &
2729 (mask_array(i,j) /= maskval)) then
2730 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2731 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2732 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2733 end if
2734 end do
2735 end do
2736 end do
2737 else
2738 do i=min_i,max_i
2739 do j=min_j,max_j
2740 do k=min_k,max_k
2741 ! Ignore masked/missing values in the source data
2742 if ((tile_array(i,j,k) /= msgval)) then
2743 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2744 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2745 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2746 end if
2747 end do
2748 end do
2749 end do
2750 end if
2752 ! end if
2754 ! Rectangle is a square of four points, and we can simply deal with each of the points
2755 else if (((max_i - min_i + 1) <= 2) .and. ((max_j - min_j + 1) <= 2)) then
2756 do i=min_i,max_i
2757 do j=min_j,max_j
2758 x_dest = where_maps_to(i,j,1)
2759 y_dest = where_maps_to(i,j,2)
2761 if (x_dest /= OUTSIDE_DOMAIN) then
2763 ! if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2764 if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2765 do k=min_k,max_k
2766 dst_array(x_dest,y_dest,k) = 0.
2767 end do
2768 end if
2770 if (present(mask_array) .and. present(mask_relational)) then
2771 do k=min_k,max_k
2772 ! Ignore masked/missing values
2773 if (tile_array(i,j,k) /= msgval) then
2774 if (mask_relational == ' ' .and. mask_array(i,j) /= maskval) then
2775 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2776 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2777 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2778 else if (mask_relational == '<' .and. mask_array(i,j) >= maskval) then
2779 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2780 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2781 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2782 else if (mask_relational == '>' .and. mask_array(i,j) <= maskval) then
2783 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2784 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2785 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2786 end if
2787 end if
2788 end do
2789 else if (present(mask_array)) then
2790 do k=min_k,max_k
2791 ! Ignore masked/missing values
2792 if ((tile_array(i,j,k) /= msgval) .and. &
2793 (mask_array(i,j) /= maskval)) then
2794 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2795 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2796 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2797 end if
2798 end do
2799 else
2800 do k=min_k,max_k
2801 ! Ignore masked/missing values
2802 if ((tile_array(i,j,k) /= msgval)) 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 do
2808 end if
2809 ! end if
2811 end if
2812 end do
2813 end do
2815 ! Not all corners map to the same grid point, and the rectangle contains more than
2816 ! four points
2817 else
2818 center_i = (max_i + min_i)/2
2819 center_j = (max_j + min_j)/2
2821 ! Recursively process lower-left rectangle
2822 call process_continuous_block(tile_array, where_maps_to, &
2823 src_min_x, src_min_y, src_min_z, &
2824 src_max_x, src_max_y, src_max_z, &
2825 min_i, min_j, min_k, &
2826 center_i, center_j, max_k, &
2827 dst_array, n, &
2828 start_x, end_x, start_y, end_y, start_z, end_z, &
2829 istagger, &
2830 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2832 if (center_i < max_i) then
2833 ! Recursively process lower-right rectangle
2834 call process_continuous_block(tile_array, where_maps_to, &
2835 src_min_x, src_min_y, src_min_z, &
2836 src_max_x, src_max_y, src_max_z, &
2837 center_i+1, min_j, min_k, max_i, &
2838 center_j, max_k, &
2839 dst_array, n, &
2840 start_x, end_x, start_y, &
2841 end_y, start_z, end_z, &
2842 istagger, &
2843 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2844 end if
2846 if (center_j < max_j) then
2847 ! Recursively process upper-left rectangle
2848 call process_continuous_block(tile_array, where_maps_to, &
2849 src_min_x, src_min_y, src_min_z, &
2850 src_max_x, src_max_y, src_max_z, &
2851 min_i, center_j+1, min_k, center_i, &
2852 max_j, max_k, &
2853 dst_array, n, &
2854 start_x, end_x, start_y, &
2855 end_y, start_z, end_z, &
2856 istagger, &
2857 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2858 end if
2860 if (center_i < max_i .and. center_j < max_j) then
2861 ! Recursively process upper-right rectangle
2862 call process_continuous_block(tile_array, where_maps_to, &
2863 src_min_x, src_min_y, src_min_z, &
2864 src_max_x, src_max_y, src_max_z, &
2865 center_i+1, center_j+1, min_k, max_i, &
2866 max_j, max_k, &
2867 dst_array, n, &
2868 start_x, end_x, start_y, &
2869 end_y, start_z, end_z, &
2870 istagger, &
2871 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2872 end if
2873 end if
2875 end subroutine process_continuous_block
2877 end module process_domain_module